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AUTOMATIZAÇÃO DO METODO DAS FORÇAS 
NA ANÁLISE DE PÓRTICOS RECTANGULARES 


J.A. TEIXEIRA DE FREITAS* e E.M.B. RIBEIRO PEREIRA ** 


SUMÁRIO 


Apresenta-se uma formulação para a análise de 
pórticos rectangulares em computador. Demonstra-se 
que o número de operações necessárias para realizar a 
análise pelo método das forças é significativamente 
interior ao exigido pelo método dos deslocamentos. 
Fornecem-se expressões analíticas para todos os opera- 
dores intervenientes e descrevem-se as principais fases 
do algoritmo. 


INTRODUÇÃO 


O método dos deslocamentos e o método das forças 
são os métodos que tradicionalmente se utilizam para 
realizar a análise elástica de estruturas reticuladas. 

O método dos deslocamentos explora o grau de 
mnideterminação cinemática da estrutura, 4. 

Bloqueando os 8 deslocamentos indeterminados q, 
determinam-se as forças Os que equilibram a solicita- 
ção. A equação do método 


K q+QG=0 (1) 


é obtida obrigando que essas forças equilibrem as pro- 
vocadas pelos deslocamentos na ausência de qualquer 
outra solicitação; K é portanto a matriz de rigidez da 
estrutura. 

O método das forças baseia-se na ideia complemen- 
tar de explorar o grau de indeterminação estática da 
estrutura, a. 


(*) Professor Aassociado, IST. 
(**) Aluno Tarefeiro, IST. 
Original recebido para publicação em 8/2/84 
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Introduzindo a libertações na estrutura, determi- 
nam-se os deslocamentos relativos que aí se desenvol- 
vem, Vo, quando se aplica a solicitação. A equação 
do método 


Fp+ty=0 (2) 


é obtida obrigando que esses deslocamentos sejam anu- 
lados pelas forças hiperstáticas correspondentes p; F é 
pertanto a matriz de flexibilidade da estrutura. 


ESCOLHA DO MÉTODO DE CÁLCULO 


É usual preterir o método das forças em favor do 
método dos deslocamentos quando se pretende realizar 
em computador a análise de estruturas reticuladas. 

Esta opção é baseada em dois argumentos. O pri- 
meiro é a facilidade de gerar automaticamaente a ma- 
triz de rigidez da estrutura, assim como o vector das 
forças nodais. O segundo é a economia que se pode 
garantir na solução numérica do sistema se se explorar 
o facto de a matriz de rigidez ser simétrica e bandeada, 
se se escolher uma numeração apropriada para os des- 
locamentos nodais. 

Como adiante se poderá verificar, é também possí- 
vel gerar automaticamente a equação de método das 
forças recorrendo a uma técnica análoga à utilizada na 
versão directa do método dos deslocamentos; a matriz 
de flexibilidade da estrutura é construída por coloca- 
ção directa da matriz de flexibilidade de cada barra 
que a compõe, e cuja definição geral se fornece ao 
programa. 

Será ainda demonstrado que, tal como a matriz de 
rigidez, a matriz de flexibilidade para além de ser simé- 
trica é também bandeada. 

Admitindo portanto que ambos os métodos são 


FIG. 1 


FIG. 2 


qualitativamente idênticos, tanto no que respeita à 
facilidade da montagem automática das equações em 
que se baseiam, como quanto à estrutura dos próprios 
sistemas, interessa finalmente comparar o tempo de cál- 
culo necessário à solução numérica do sistema resolu- 
tivo associado a cada um destes dois métodos. 

Como termo de comparação, considerem-se pórti- 
cos rectangulares do tipo do representado na figura 1. 

Se V representar o número de vãos, P; o número 
de pisos do vão k e 


Pi= mix (P,s Pa) 


conclui-se que 


V-l V 


B=PL+H+HPr+ZP+ a, 
k=1 ki 
VA 
N=P+P+2ZP, 
ke 


em que B é o número de barras e N o número de nós 
livres. 

O grau de indeterminação cinemática da estrutura 
é definido por 


B=3N 
e o grau de indeterminação estática por 
q =3(B-—N) 


As expressões acima mostram que a dimensão do 
sistema do método das forças é sempre inferior 
à do método dos deslocamentos (B < 2 N). 

A dimensão da semi-banda é definida por 


b<3+3.min (max (Py),V) 


em que VC = V para a equação do método das forças 
cV'=V+1 paraa do método dos deslocamentos. 
A dimensão da semi-banda da matriz de-flexibili- 
dade nunca é, portanto, superior à dimensão da semi- 
banda da matriz de rigidez da mesma estrutura. 
Como o número de operações numéricas necessá- 


rias à solução de um sistema simétrico de dimensão D e 
semi-banda b é definido |1| por 


T=(3D-2b) b”/6 (3) 


concluiu-se que, para um dado pórtico, rectangular, o 
tempo dispendido na solução da equação do método 
dos deslocamentos é sempre superior ao necessário 
para a solução da equação do método das forças. Esta 
diferença é ainda maior se a rotina utilizada não explo- 
rar o bandeamento da matriz, situação em que na 
equação (3) a dimensão da maior semi-banda iguala a 
dimensão do sistema (b = D). 

No caso particular dos pórticos regulares como o 
representado na figura 2 (P' = P, sendo P o número 
de pisos) as definições acima formam as seguintes 
expressões 


€ = 3PY,Ê=3PIVA 1) 


A relação entre o número de operações necessárias 
para a solução das equações do método dos desloca- 
mentos, T 9º & do método das forças, T , está repre- 


sentada nos diagramas das figs. 3 e 4 para vários valo- 
res do número de vãos e de pisos em que R=To/T ; 
a 
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EQUAÇÃO DO METODO DAS FORÇAS 
Seja 
= Mm = 1,2, (4) 


o vector dos esforços independentes, isto é, os esforços 
necessários e suficientes para caracterizar o estado de 
tensão em qualquer secção da estrutura. 

Tornando a estrutura isostática através da introdu- 
ção de a libertações onde se aplicam as forças hipers- 
táticas 


p = (pr), dk = (5) 


encontra-se a seguinte definição para condição de equi- 
líbrio da barra m 

> 5 — Bm p + Aam m =— Eai B (6) 
em que a coluna j da matriz B, representa a distribui- 
ção de esforços que equilibra a força hiperstática p; = 
= | e Xom é O vector que lista os esforços causados 
pelo carregamento. 


Seja 
V = (vi 5 VD genes Vs (7) 
o vector dos deslocamentos correspondentes às forças 
hiperstáticas e 


U = (My Usos dp) (8) 


o vector das deformações independentes corresponden- 
tes aos esforços (4). 

Na hipótese das pequenas deformações e desloca- 
mentos, o deslocamento v; introduzido na estrutura 
hiperstática pela acção das deformações u, é definida 
por 


| 


MM 


v= 


Ba Um (9) 


ma. 


h 


em que Bm representa a transposta da matriz Ba. 

Se se admitir que as peças estruturais têm um com- 
portamento elástico-linear, as relações que caracteri- 
zam as deformações no elemento m podem ser expres- 
sas na forma 


u=FX.+u == 1d B (10) 


em que F' representa a matriz de flexibilidade do ele- 
mento é Um O vector que lista as deformações causa- 
das pela solicitação de vão. 

Agrupando (6), (9) e (10) de acordo com as defi- 
nicões (4), (5), (7) e (8) encontram-se as seguintes 
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expressões para as condições de equilíbrio, compati- 
bilidade e elasticidade da estrutura isostática 


Xx=Bp+X,v=Bu,u=F X+uw (11-13) 


em que P=[Folsti= 1,220 B 


é uma matriz diagonal por blocos e 


B -—— | Bm | + mMm= l 2 PSI AS B 

Para tornar as condições (11-13) válidas para a 
estrutura hiperstática a analisar, basta impor que sejam 
nulos os deslocamentos relativos onde se introduziram 
as a libertações para obter a estrutura isostática: 


v=0 (14) 


A equação (2) do método das forças, em que 


F=BFB,w=B(FX +u) (15,16) 


é obtida substituindo as relações de elasticidade (13) 
na condição de compatibilidade (12, 14) e eliminando 
os esforços na equação resultante através da condição 
de equilíbrio (11). 

Em consequência da matriz de flexibilidade das 
barras ser diagonal por blocos, os definições (15, 16) 
podem ser escritas na forma 


B B 
F= X Fm,VYo= E VYom (17, 18) 
em que Fa = Bo E. Ba (19) 
= Yom =— Bm Lim (20) 
com Uma EF Xa É Um (21) 


DEFINIÇÃO DOS OPERADORES DE CAUSALI- 
DADE ELÁSTICA 


Como se ilustra na figura 5, os esforços escolhidos 
para caracterizar o estado de tensão no elemento m 


im — (M, , Ms , NJm (22) 
são os momentos flectores nas secções extremas e o 
estorço axial na extremidade da direita, para a orien- 


tação ai indicada. 


Analogamente, o vector das deformações generali- 
zadas 


Um = (01,0>,€)m (23) 


agrupa as rotações das secções extremas, medidas em 
relação à corda do elemento, e o encurtamento axial 
por ele sofrido, como se indica na figura 6. 

Se se admitir que a barra m é homogénea, de eixo 
recto e de secção constante, com momento de inércia 
Im € área Am, encontra-se a seguinte expressão para a 
matriz de flexibilidade 


F,=|L/3EL L/6EL O 
' L/3EL O (24) 
simétrica L/EA 


a qm 
que define as deformações generalizadas (23) causadas 
pelos esforços correspondentes (22); na definição (24) 
L representa o comprimento inicial da peça e Em O 
seu módulo de elasticidade. 


(Fig. 6) 


As deformações adicionais são definidas por 


L? 
adafo x 
a a). 


ab +b (1-3b)/L] | M | 


DR RR 


— 6al/AL) À; 


cas E la mo HT 
para as forças concentradas apresentadas na figura 5. 


Outras acções de vão, nomeadamente forças distri- 
buídas, variações de temperatura e pré-esforço, são con- 
sideradas em [2], podendo encontrar-se na mesma refe- 
rência a generalização das definições (24) e (25) para 
peças de secção variável. 


DEFINIÇÃO DOS OPERADORES DE EQUILÍBRIO 


Na condição de equilíbrio (11), X, representa uma 
qualquer distribuição de esforços que equilibra o carre- 
gamento, globalmente ou cada uma das acções tomadas 
separadamente [3]. 

Analogamente, as colunas da matriz B representam 
as distribuições de esforços que equilibram as forças 
hiperstáticas correspondentes; por definição são a dis- 
tiibuições de esforços auto-equilibradas e linearmente 
independentes [4]. 

É em conseguência destas propriedades que se 
podem estabelecer as expressões gerais para os opera- 
dores de equilíbrio que a seguir se apresentam. 


Acção das Forças Hiperstáticas 

Como se ilustra nas figs. 7 e 8, um pórtico rectan- 
gular pode ser considerado como tendo sido construído 
ligando entre si 


C=B-N (26) 


quadros rectangulares; o grau de indeterminação está- 
tica da estrutura é portanto 


“= 35€ (27) 


(Fig. 7) (Fig. 8) 


Considere-se o quadro genérico M representado na 
figura 9, onde se indica a numeração utilizada para 
identificar as barras que a constituem, assim como a 
orientação adoptada. 

Para tornar o quadro isostático introduz-se um 
corte, escolhendo-se para forças hiperstáticas 


Pu = (pi,p>,Psm (28) 


um sistema ortogonal actuando no centro geométrico 
do quadro, ligado ao corte introduzido por meio de 
bielas rígidas, como se ilustra na figura 10. 
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" 2 


(Fig. 9) 


(Fig. 10) 


As forças hiperstáticas provocam uma distribuição 
de esforços na barra i definida por 


Xim = Bim pu,i=1,2,3,4 (29) 
em que 
Be=| 1-1 i[Bas] t=t =i7 
af cad 11 1 
| 2/1 ) E 2/Li 
Ba -t 2 =1 [hu=f=t at 47 
si =] e =. sl | 
—-2/L|. ova . 
dovy (30) 


se se fizer pp = 1,p: = 2/L, e P; = 2/L;; nas defini- 
ções acima Lm representa o vão da barra m. 


Acção do Carregamento 


O processo mais simples para equilibrar forças apli- 
cadas a um nó da estrutura consiste em descarregá-las 
na consola que se obtém desligando do resto da estru- 
tura o montante que termina nesse nó, como se ilustra 
na figura 11. 


Os esforços no pilar m desse montante são defini- 
dos por 


Mom — Bomn a (31) 
com Bom =[|y, O +1| 
y 0 +1 (324) 
0 1 à 
em que W=)—Yo k=i,j 


devendo notar-se que 
Dom = 0 (32b) 


se y; for positivo. 
As forças aplicadas nas barras da estrutura são 


equilibradas transformando numa consola a malha a 
que essas barras pertencem, como se ilustra na figura 12. 
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Os esforços nas barras que pertencem à malha são 
definidos por 


Xm= UU Mb AM —h A; — AM; AD (334) 


Xom = = 4,20, 0, (33b) 
XKom = Xom = O (33c,d) 


se se adoptar a numeração e a orientação indicada 
na figura 9, 


As reacções de apoio são definidas por 


= MAM MAMMA YA) 


Asm 


A2m 
(Fig. 12) 


e introduzem esforços definidos por (31) nos pilares 
pertencentes a esse montante. 

Nas definições acima, X/ pode representar a resul- 
tante de uma força distribuída. 


DESCRIÇÃO GERAL DO ALGORITMO 


Tendo definido as expressões gerais dos operado- 
res intervenientes na análise de pórticos rectangulares 
pelo método das forças, interessa descrever sumaria- 
mente as principais fases do algoritmo de solução a 
adoptar na preparação de um programa de cálculo 
automático |5|. 


Caracterização da Estrutura 


A regularidade da geometria dos pórticos rectan- 
gulares permite minimizar a intervenção do utilizador 
no programa de cálculo automático. 

Como dados o utilizador tem de fornecer ao pro- 
grama o número de vãos, o número de pisos por vão, 
a largura de cada vão e a altura de cada piso. Com 
base nesta informação, podem construir-se rotinas que 
numeram automaticamente os nós, os quadros e as 
barras da estrutura e que definem para cada barra os 
nós que a limitam e os quadros a que pertence. A nu- 
meração dos quadros deve ser escolhida de modo a 
garantir o bandeamento da matriz de flexibilidade. 

Para caracterizar as propriedades mecânicas da 
estrutura, o utilizador tem ainda de fornecer a área, O 
momento de inércia e o módulo de elasticidade para 
cada grupo de barras que partilhem as mesmas carac- 
terísticas. 


Caracterização do Carregamento 


O último grupo de dados que o utilizador tem de 
fornecer é o que caracteriza a solicitação, indicando o 
tipo, a localização, a intensidade e o sentido. 

Definidas as forças nodais determinam-se as entra- 
das correspondentes do vector dos esforços iniciais, Xo, 
recorrendo à definição (31). 

Com base na informação sobre a solicitação de vão, 
calcula-se a sua contribuição para os esforços iniciais 
recorrendo a (33) e (31) e as deformações iniciais que 
provoca utilizando a definição (25). 


Montagem do Sistema Resolutivo 
Seja Bm = [Bm], n=1,2,..., a 


a expressão geral da matriz de equilíbrio presente nas 
definições (19) e (20). 

Como uma barra pertence no máximo a dois qua- 
dros, pelo menos « — 2 e no máximo « — 1 das sub- 
matrizes Bmn são nulas; o bandeamento da matriz de 


flexibilidade (17) é uma consequência desta proprie- 
dade. 

O termo genérico da matriz de flexibilidade toma 
expressão 


Fy = À, Bmi Fm Bm (34) 


estendendo-se o somatório às barras m que pertencem 
simultaneamente aos quadros i e j, encontrando-se ana- 
logamente a seguinte definição para o termo genérico 
do vector dos deslocamentos nas libertações 


B 7 Qua q + Un) (35) 


à = 3 
mEi 


Nas expressões (34) e (35) os operadores de equilí- 
brio e de elasticidade associados à barra m são defini- 
dos por (30) e (24), respectivamente. 

Depois de montado o sistema resolutivo (2), levan- 
ta-se a hiperstaticidade da estrutura recorrendo a roti- 
nas que explorem a simetria e o bandeamento da matriz 
de flexibilidade. 


ESTRUTURAS COM APARELHOS DE LIBERTA- 
ÇÃO 


A formulação anteriormente apresentada foi desen- 
volvida supondo que a estrutura a analisar era consti- 
tuída por barras que se ligavam rigidamente entre si e 
ao meio de fundação. Indicam-se em seguida as altera- 
rações que se devem introduzir nessa formulação 
quando algumas dessas ligações se realizam por meio 
de aparelhos de libertação como os representados na 
figura 13. 


Libertações Elásticas 


Seja X =X) m= 12, B 


o vector que agrupa as forças que se desenvolvem nos 
B" aparelhos de libertação e 


+ ? 7 
=>) ti=l2a B 
A lis F 
com e do (36) 
o vector dos deslocamentos correspondentes; 


PE m= 2 B 


TT 


61 


é a matriz diagonal que agrupa as constantes elásticas 
das molas. 

Conhecendo o tipo e a localização dos aparelhos 
de libertação é fácil exprimir as forças nas molas em 
função dos esforços independentes das barras em que 
esses aparelhos se situam; a condição de equilíbrio 


X =B.p+Xo (57) 


pode portanto ser obtida por particularização da defi- 
nição (6) para os esforços independentes. 

Atendendo a que a condição de compatibilidade 
(9) toma agora a forma 


B Ex 
v=2Bu + 4XB ul (38) 
mm] mãI 


encontra-se a seguinte generalização 


F=BEFBA+B Fo B' (39a) 


,=B(FX+u)+B FX, (59) 


para as definições (15, 16) dos operadores presentes 
na equação do método das forças. Os termos adicionais 
podem ser obtidos usando as técnicas anteriormente 
referidas, pois os operadores intervenientes mantém as 
mesmas propriedades, 


Libertações Perfeitas 


Para as B” libertações livres que existam na estru- 
tura podem ser estabelecidas equações de equilíbrio do 
tipc da das libertações elásticas, 


X=0=Bºp+X” 


m Hm 


m= 12. 56 440) 


devendo agora impor-se que estas forças sejam nulas 
em substituição das relações de elasticidade (36). 

Como a condição de compatibilidade (38) se gene- 
raliza agora para 


B DM. Bº. 
v=XBu+2ZB uv +XB/uí 
mel mal ml 


encontra-se a seguinte expressão para o sistema do 
método das forças 


F|BIlpl+| vw | =o 


B” ur we” 


(41b) 


em que a segunda equação representa a condição de 
equilíbrio (40). As definições (39) são ainda válidas 
para a equação da compatibilidade (41a). 
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Como o sistema (41) mostra, o grau de hiperstati- 
cidade da estrutura passa a tomar o valor 


[o Má = -— BY 


como aliás seria de esperar. 

O método acima descrito deve ser abandonado 
quando o número de libertações perfeitas se torna 
significativo relativamente ao grau de hiperstaticidade 
da estrutura na ausência dessas libertações. Nestes 
casos está-se a resolver uma estrutura de hiperstatici- 
dade a” <a utilizando um sistema resolutivo de dimen- 
são «a + Bº>a. 

Um método alternativo consiste em modificar direc- 
tamente a matriz de flexibilidade das barras para se 
passar a contemplar o efeito das libertações. 


Seja XM = By Pm + XoM (42) 


a equação de equilíbrio dos quadros que partilham a 
barra onde existem libertações perfeitas; o sistema (42) 
é obtido agrupando as condições (6) para os quadros 
em questão. 

A equação (42) pode ser expressa na forma equiva- 
lente 


| Xo | (434) 
x” | (43b) 
) Es + M 


B. | | p' + 


x =olilb |b lp 


— — 0 mu 


em que se explicitaram os esforços que devem ser nulos. 
O subvector p, é escolhido de modo a garantir que a 
submatriz b;m seja não-singular. 

Resolvendo a equação (43b), encontra-se a defini- 
ção alternativa 


Xu pre By pa É AM 
para a condição de equilibrio (42), em que 
By - (B, — B, b; b, Dm 


C x =(X— Bb; Ao dus 


ol ú 


RESOLUÇÃO DA ESTRUTURA 


Levantada a hiperstaticidade da estrutura, os esfor- 
cvs Independentes e as forças nas libertações elásticas 
são calculadas a partir das condições de equilíbrio (6) 
e (57), respectivamente; de entre as forças hiperstá- 
ticas, só contribuem as correspondes aos quadros que 
partilham a barra em questão. 

Os esforços dependentes podem depois ser calcula- 
dos a partir das seguintes definições 


Nº=N-M,T=T+» 
Tr =(—M+M-—aLh+hM)/L 


de acordo com a notação indicada na figura 5. 

De acordo com a condição de equilíbrio (31), os 
deslocamentos correspondentes às forças nodais repre- 
sentadas na figura 11, são definidas por 


= 2 


mfPn 


| Bo - Vad 


(da + BI um ] 
devendo o somatório abranger apenas os pilares P, que 
ligam o nó n à fundação. As deformações independen- 
tes, u, e os deslocamentos nas libertações elásticas, u”, 
são calculadas a partir das relações constitutivas (10) e 
(36), respectivamente, sendo os deslocamentos nas 
libertações perfeitas fornecidos pela solução do siste- 
ma (41). 


CONCLUSÃO 


Tal como sucede na equação do método dos deslo- 
camentos, a matriz do sistema do método das forças 
é simétrica e bandeada. Em qualquer dos casos a mon- 
tagem automática dessa matriz é facilmente realizável, 
apresentando o método das forças a vantagem de pro- 
duzir sistemas cuja solução envolve um menor volume 
de operações numéricas. 

A viabilidade de implementar o método das forças 


em computador depende essencialmente da possibili- 
dade de gerar automaticamente uma distribuição de 
esforços que equilibre a solicitação. Para o caso parti- 
cular dos pórticos rectangulares esta operação pode ser 
facilmente executada explorando a regularidade da 
topografia da estrutura. 

A existência de libertações perfeitas complica tanto 
a montagem como a solução do sistema resolutivo, 
qualquer que seja o método de análise. Tal como no 
método dos deslocamentos, também no método das for- 
ças se pode recorrer a uma de duas técnicas alternativas 
para processar o efeito das libertações perfeitas; aceitar 
os deslocamentos nessas libertações como variáveis do 
sistema ,ou alterar directamente as condições de equili- 
brio e de compatibilidade das barras afectadas por 
essas libertações. A primeira destas técnicas tem a van- 
tagem de exigir um esforço mínimo de programação 
adicional, 
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SIMULAÇÃO COMPUTACIONAL DO GOLPE DE 
ARÍETE EM SISTEMAS ELEVATÓRIOS COM 
CONDUTAS DE ASPIRAÇÃO LONGAS 


A. BETÂMIO DE ALMEIDA* 
P. FALCÃO E CUNHA** 


RESUMO 


Em alguns casos de sistemas elevatórios a conduta de 
aspiração tem um comprimento significativo e os efeitos do 
golpe de aríete podem ser alterados de forma apreciável 
relativamente aos casos usuais de sistemas com grupos elec- 
trobomba colocados junto do reservatório de montante. 

No presente artigo faz-se a síntese dos principais resulta- 
dos obtidos com a exploração sistemática de um modelo 
matemático desenvolvido para o efeito, apresentando-se 
algumas conclusões respeitantes a cotas piezométricas ex- 
tremas e respectivas envolventes, ao tempo de anulação de 
caudal nos grupos e à velocidade residual de rotação daque- 
les. Estas conclusões são apresentadas de forma adimen- 
sional com vista à sua utilização prática e à generalização 
da sua aplicabilidade. 


ABSTRACT 


In some pumping systems the suction line is significantly 
long and the water hammer effects can be changed in some 
extent in relation to the usual situations where the pumps 
are placed near the upstream reservoir. 

In this article a synthesis is made of the main results 
obtained with the systematic operation of a mathematical 
model develloped to cope with such situations. The authors 
present the conclusions concerning extreme pressure heads 
and respective envelope lines, time of flow stoppage in the 
pumps and residual pump speed. 

These conclusions are presented in an adimensional wav 
aiming its practical application and the generalization of its 
field of utilization. 


PARTE | — Modelação Matemática 


1. INTRODUÇÃO 


Nos casos usuais de sistemas elevatórios, os grupos elec- 
trobomba estão colocados junto ao reservatório de mon- 
tante sendo o comprimento das respectivas condutas de 
aspiração, para efeitos de análise do golpe de aríete, des- 
prezável em comparação com o comprimento da conduta 
elevatória a jusante dos grupos. Para esta situação existem 
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já numerosos elementos de consulta que possibilitam a de- 
terminação aproximada das cotas piezométricas extremas 
que ocorrem após a saida de serviço simultânea dos grupos 
electrobomba nos casos de condutas únicas e uniformes, 
sem dispositivos especiais de protecção contra o golpe de 
ariete e contando-se unicamente com a Inércia dos grupos e 
de eventuais massas adicionais acopladas (volantes de inér- 
cia). Referências a estes elementos podem ser encontrados 
em ALMEIDA, 1981 e 1982. 

Em alguns casos, contudo, a conduta de aspiração tem 
um comprimento significativo e os efeitos do golpe de 
ariete podem ser alterados de forma apreciável. Este facto 
pode eventualmente motivar o projectista a considerar dife- 
rentes localizações da central elevatória ao longo da con- 
duta, tendo em conta as características topográficas e 
outros condicionalismos, por forma a conseguir atenuar as 
variações extremas de pressão no sistema e obter uma solu- 
ção globalmente mais econômica. 

Na parte | do artigo apresentam-se os aspectos principais 
da modelação matemática e do programa de cálculo auto- 
mático (denominado RESBOM) desenvolvido, bem como 
uma breve análise qualitativa do fenómeno. Na parte II do 
artigo apresenta-se uma sintese dos resultados obtidos com 
a exploração do programa RESBOM para simulações de 
regimes transitórios decorrentes da saida de serviço de gru- 
pos electrobomba colocados, em paralelo, em diversas sec- 
ções ao longo de uma conduta elevatória uniforme. Proce- 
de-se, também, a uma análise geral do problema baseada 
em parâmetros adimensionais cuja dedução se apresenta na 
parte II, Os resultados obtidos incluem as envolventes das 
cotas piezométricas extremas, o tempo de anulação do 
caudal nos grupos a velocidade de rotação dos grupos no 
instante em que se dá a anulação do caudal. Uma breve 
análise de sensibilidade, relativa à variação do passo de 
cálculo À t e à perda de carga hidráulica, é também apre- 
sentada na parte II do artigo. 

É de referir ser escassa a informação bibliográfica acerca 
deste assunto, sendo de salientar como de maior interesse 
prático uma comunicação de Lev Pavluch (PAVLUCH, 
1970). 


* Eng. Civil, Prof. Auxiliar do IST, Membro do CEHIDRO. 
** Eng. Civil, Ex. Assistente do IST, Membro do CEHIDRO. 


Original recebido para publicação em 19/6/84. 


-—— — LINHA PIEZOMETRICA ESTÁTICA 

—— - LINHA PIEZOMEÉTRICA EM REGIME PERMAMENTE 
— — ENVOLVENTE PIEZOMÉTRICA MÁXIMA 

—— ENVOLVENTE PIEZOMÉTRICA MÍNIMA 


FIGURA 1 — ESQUEMA GERAL DE UM SISTEMA ELEVATÓRIO 


2. ANÁLISE QUALITATIVA 


Após a saída de serviço simultânea dos grupos electro- 
bomba (situação clássica para análise do golpe de ariete), 
desencadeia-se um regime transitório que difere daquele 
que ocorreria se os mesmos grupos estivessem colocados 
junto do reservatório de montante (Figura 1). De modo 
qualitativo, as características principais deste regime são as 
seguintes: 

— no trecho da conduta a montante dos grupos (con- 
duta de aspiração) a linha piezométrica começará por 
subir, em consequência da diminuição do caudal na 
respectiva extremidade de jusante (bombas); 

— no trecho de conduta a jusante dos grupos a linha 
prezométrica começará por baixar, em consequência 
da diminuição do caudal na respectiva extremidade 
de montante (bombas), tendendo as variações extre- 
mas de pressão neste trecho a serem inferiores às que 
ocorreriam no caso de a conduta de aspiração ter um 
comprimento nulo; 


— a variação das cotas piezométricas instantâneas nas 
secções de entrada e de saida das bombas poderá pro- 
vocar a passagem do escoamento através destas e al- 
terar o instante de anulação do caudal ou de fecha- 
mento das válvulas de retenção. 


Após o fechamento das válvulas de retenção colocadas a 
jusante dos grupos, os regimes transitórios nos dois trechos 
da conduta elevatória ficam desacopolados, oscilando as 
linhas piezométricas respectivas em torno da posição de 
equilíbrio final. Na figura 2 exemplifica-se a variação das 
cotas piezométricas nos primeiros instantes após a saída de 
serviço dos grupos. 


3. HIPÓTESES E EQUAÇÕES BÁSICAS 
3.1 — Hipóteses básicas 


Na elaboração do programa RESBOM admitiram-se as 
seguintes hipóteses: 


— conduta elevatória única e uniforme (diâmetro, espes- 
sura, material e celeridade constantes ao longo de 
todo o comprimento); 

— saída de serviço simultânea dos grupos electrobomba 
podendo os mesmos ser substituídos por um grupo 
equivalente; 

— ausência de dispositivos especiais de protecção à excep- 
ção das válvulas de retenção colocadas imediata- 
mente a jusante dos grupos e de eventuais volantes de 
inércia acoplados a estes. 

— fechamento instantâneo das válvulas de retenção nos 
instante da anulação do caudal nas bombas, admitin- 
do-se que as mesmas não reabrem. 


Na dedução das equações básicas do modelo matemático 
admitiram-se ainda as seguintes hipóteses: 


e —— 
mi | 


FIGURA 2 = SISTEMA ELEVATÓRIO COM A BOMBA NUMA SECÇÃO IMTERMÉDIA EVOLUÇÃO DA 


LINHA PIEZOMÉTRICA APÓS SAÍDA DE SERVIÇO DOS GRUPOS(sem perda de carga! 
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— escoamento puramente turbulento e unidimensional 
nas condutas, sendo desprezáveis a altura cinética e as 
perdas de carga localizadas e admitindo-se, em cada 
secção da conduta, a distribuição uniforme das velo- 
cidades; 

— validade, durante os regimes transitórios, das curvas 
características e das fórmulas das perdas de carga 
contínuas utilizadas em regime permanente (hipótese 
quasi-estacionária); 

— comportamento reológico elástico e lincar do liquido 
(água) e da conduta, não sendo considerados os even- 
tuais efeitos de cavitação ou de rotura da coluna li- 
quida. 


3.2 — Equações básicas 

As equações básicas, válidas na modelação dos transitó- 
rios hidráulicos ao longo da conduta, são as seguintes 
(CHAUDHRY, 1979; WYLIE e STREETER, 1978): 

— Equação do Balanço da Quantidade de Movimento 


dQ 
Õ! 


+ RQIQ| = 0 (1) 


9H 
"Ee Dx 


— Equação da Continuidade 


OH o é dO 
+ ES as = À (2) 


e] 
+ 


em que Q = caudal, H = cota piezométrica, S = área da 
secção transversal da conduta, R = coeficiente de perda de 
carga unitária, c = celeridade das ondas elásticas, g = ace- 
leração da gravidade, x = distância ao longo do eixo da 
conduta et = tempo. 

Na resolução numérica do sistema de equações (1) e (2) 
utiliza-se o método das características, na modalidade de 
malha computacional prefixada com À t constante, em 
conjunto com o método explícito de diferenças finitas. Da 
aplicação deste método misto resultam as seguintes duas 
familias de equações (ALMEIDA, 1981): 


— Equações Es 
| c 
tHpo— Hg) + gs (Op — MB) + 


+ RQg |Qgl (xp — xp)=0 (3) 


Xp — Xxp=+cÃt (4) 
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— Equações CC 
| c 
Ho — Hn) -— (Qp — 
(Hp D) ES (Qp — Qp) + 


+ RQp |Opl xp — xp)=0 (5) 


em que os indices indicam os pontos da malha computacional 
no plano (x.t) a que dizem respeito os valores (Figura 3). 
Nocasode 4x = c At, os pontos E e D coincidem com 
os nós A e B da malha prefixada, respectivamente à es- 
querda e à direita do ponto de cálculo P e no instante de 
cálculo anterior t — À t. No caso de ser escolhido um 
passo de cálculo temporal tal que A x > c&át, os va- 
lores de Q e de H correspondentes aos pontos E e D são 
obtidos por interpolação linear a partir dos valores nos nós 
A, Be € (Figura 3). A execução de interpolações lineares 
dá origem a um amortecimento de tipo numérico cujo 
efeito há que limitar e deve estar presente na escolha dos 
valores de A x e de A t. A resolução das equações 
(3) e (5) permite obter de forma explícita os valores de Qp 
e de Hp nas secções internas das condutas e em cada ins- 
tante de cálculo. 


3.3 — Modelação dos grupos electrobomba 


A obtenção de soluções numéricas do sistema de equa- 
ções básicas exige a caracterização matemática das condi- 
ções de fronteira e o conhecimento das condições iniciais. 

No caso em análise houve que considerar dois tipos de 
condições de fronteira (Figura 3): grupo electrobomba 
equivalente (na extremidade de montante ou numa secção 
intermédia da conduta elevatória) e reservatórios de gran- 
des dimensões (nas extremidades de montante e de jusante 
da conduta elevatória). 


Em regime permanente, a caracterização do grupo elec- 
trobomba equivalente a n bombas iguais e colocadas em 
paralelo é feita unicamente com a equação da curva carac- 
terística das bombas a qual se admite ser, no programa 
RESBOM, do seguinte tipo: 


Hp = ANº + BNQp /n — CQh ?/n? (7) 


= E 
E 3 

“| Ed 

| = = 
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ida at 
= [a 


FIGURA 3 — MODELAÇÃO MATEMÁTICA (método das carmeteristicas) 
DISCRETIZAÇÃO E CONDIÇÕES DE FRONTEIRA 


em que Hy = altura de elevação nas bombas (diferença 
entre as cotas piezométricas imediatamente a montante e a 
jusante das bombas), Qp = caudal no grupo equivalente, 
N = velocidade do grupo equivalente (em r.p.m.), A, B, 
C = coeficientes da curva característica de cada bomba 
fornecida pelo fabricante. A equação pode ser adaptada de 
modo a caracterizar bombas colocadas em paralelo mas 
diferentes entre si ou bombas colocadas imediatamente em 
série (FOX, 1977). 

Desligados os motores eléctricos e iniciado o regime 
transitório, o valor de N diminuirá. No programa RES- 
BOM, a determinação do valor de N em cada instante de 
cálculo é feita de dois modos diferentes: 


— No intervalo de tempo Tp (tempo de reflexão elás- 
tica total do sistema) após a saída de serviço dos gru- 
pos electrobomba segue-se um modelo conceptual de 
tipo exponencial baseado numa solução exacta apro- 
ximada (ALMEIDA, 1981). 

— Para tempos de simulação do regime transitório su- 
periores a Tp procede-se à integração numérica da 
equação das massas girantes de acordo com uma téc- 
nica de previsão-correcção. 


As equações utilizadas na modelação matemática são as 
seguintes: 


— Para À T& Te 


B À 


Me Meda ste E ya (8) 


com B=1/(mn,A Te) 


— Para À T > Tg 
N=N ar -05KgE(M+M, Ap) dt (9) 
com Kg = 3600g / (HI? PD?) 


sendo À t = Intervalo de tempo após o início do regime 
transitório, Tg = tempo de reflexão total (Tp = 2L,/c), 
No N, ar = valoresde N nos instantes de cálculo tet 
— À t, n, = rendimento do grupo electrobomba, p = pa- 
râmetro característico de conduta (ver Anexo |), À = coe- 
ficiente de compensação energética (ver Anexo 1), PD? = 
parâmetro característico da inércia das massas girantes e 
M, M, a, = binário actuante no grupo electrobomba 
calculado do seguinte modo, sendo y = peso volúmico da 
água estando incluída em Kg a conversão numérica de uni- 
dades: 


M =y Qb Hb / (mo N) (10) 


Na fase de exploração do programa a que o presente 
estudo diz respeito, admitiu-se que os coeficientes A, Be C 
da curva característica da bomba permanecem constantes 
durante o regime transitório e que a utilização das equações 
(9) e (10) é válida até ao instante de fechamento da vál- 
vula de retenção do grupo electrobomba. Não é assim feita 
a distinção entre as diferentes fases de funcionamento pos- 


síveis da bomba: fase de funcionamento normal, fase de 
dissipação e fase de turbinagem. São impostas, contudo, 
duas condições adicionais correspondentes às duas últimas 
fases: 


— a velocidade de rotação não pode ser negativa (valor 
mínimo = Ir.p.m.); 

— em cada instante de cálculo, a velocidade de rotação 
N, não pode ser superiora 10 N, as 


Uma caracterização mais correcta do comportamento 
das bombas hidráulicas exigiria o conhecimento das curvas 
características (alturas de elevação e binários) nas eventuais 
fases de funcionamento em turbinagem e em dissipação o 
que é, em geral, muito dificil de obter junto dos fornecedo- 
res. 

Da conjugação das equações do método das caracteristi- 
cas (equações (3) e (5) ) com as equações (7) e (8) ou 
(9) obtém-se a seguinte expressão para o caudal instantá- 
neo na bomba equivalente: 


de= tyr xy y =4x0)/2]/€ (11) 


sendo: 
y=BN,/n —- a c/(gS) 
k= Cl — C2 — AN (12) 
É =C/ nº 


em que (Figura 1), no caso de bombas de extremidade 
(L, = 0), a = le Cl = Zy e no caso de bombas com 
conduta de aspiração longa (L| * 0), a =2e 


Ci = Hp + (— 


s — RIQrl). Qr (13) 


sendo, em ambos os casos, 
“RIQ5I).Qp (14) 


Com base na expressão (11) o programa detecta automa- 
ticamente a tendência de inversão de sentido do caudal Qp 
e simula o fechamento da válvula de retenção (Q, = 0). 

Determinado o valor de Qp em cada instante de cálculo, 
os valores das cotas piezométricas imediatamente a montante 
e a jusante do grupo electrobomba equivalente são obtidos 
a partir das equações (3) e (5). 


3.4 — Reservatórios de extremidade 


Admite-se, no presente estudo, a existência de reservató- 
rios de grandes dimensões (Figura 3) nos quais a cota pie- 
zométrica coincide com a cota da respectiva superfície livre. 
Nos casos em que a conduta de aspiração tem um compri- 
mento não nulo (L, * 0), a condição de fronteira na 
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secção inicial da conduta (extremidade de montante) é mo- 
delada substituindo Hp. na equação (5), por Zy é obten- 
do-se assim directamente a partir desta equação o corres- 
pondente valor de Qp. No caso de a conduta de aspiração 
ter comprimento nulo (L| = 0) a modelação do reservató- 
rio de montante está implícita na modelação dos grupos 
electrobomba (expressões (11) e (12) ). 

Na extremidade de jusante, a condição de fronteira é 
modelada de modo análogo substituindo Hp, na equação 
(3), por Z, o que permite o cálculo directo de Qp. 


4. PASSO DE CÁLCULO TEMPORAL 


Com base nos valores fornecidos de L,, Ls e do número 
de trechos de cálculo NT, o programa RESBOM determina 
automaticamente o valor do intervalo de tempo At * com- 
patível com a condição de estabilidade numérica: 


L, + L» 
rd = 2 (15) 


Para diminuir os erros na integração da equação das 
massas girantes (equação (9) ), o programa calcula o 
passo de cálculo A t”” definido por (ALMEIDA, 1981): 


At" =DTpAÃ/2 (16) 


e adopta para passo de cálculo temporal o menor dos dois 
valores de À t. À expressão (16) corresponde à condição 
de a variação de velocidade do grupo electrobomba não ser 
superior, num esquema de integração explicito de 1º or- 
dem, ao valor inicial de N. 

O tempo total de cada simulação é automaticamente es- 
colhido por forma a que a envolventes possam incluir os 
valores das cotas piezométricas máximas que ocorrem após 
o fechamento das válvulas de retenção. Admitindo que a 
simulação de regime transitório se inicia no instante t = 0, 
o instante final da simulação, Tp, é calculado a partir da 
seguintes expressão (ALMEIDA, 1981): 


pe) vs 


TeE()=Tw(1+A)+20 (17) 


em que Tw = tempo de inércia rigida da coluna liquida 
(Tw = 0,5 p Tg) 


5. EFEITOS DA TÉCNICA DE INTERPOLAÇÃO 


A adopção de um passo de cálculo A t inferior ao valor 
à 1º, calculado pela expressão (15), exige a execução de 
interpolações lineares ao longo do eixo da conduta no cál- 
culo de Hp e de Op o que introduz alterações ou erros na 
simulação numérica que são do seguinte tipo (ALMEIDA, 
1981): 


— alteração nas amplitudes das flutuações de pressão 


(amortecimento numérico fictício); 


— alteração na velocidade de propagação ou de fase. 


68 


A análise completa destes efeitos é muito dificil, em espe- 
cial no caso de ocorrerem múltiplas reflexões das ondas 
elásticas e quando as perdas de carga hidráulica não são 
desprezáveis. Em casos muito simplificados é possível a 
aplicação da técnica de análise de Fourier, nos domínios do 
espaço e do tempo, conseguindo-se então obter os factores 
de atenuação e de propagação de ondas sinusoidais em 
função dos respectivos comprimentos de onda (WIGGERT 
e SUNDQUIST, 1977). Do ponto de vista prático, os efei- 
tos motivados pela técnica de interpolação podem ser 
muito importantes pois podem conduzir a cotas piezomé- 
tricas máximas calculadas inferiores às reais e a velocidades 
ficticias de propagação das perturbações, ao longo das 
condutas, superiores às reais. 

Na simulação de regimes transitórios numa conduta ele- 
vatória, a importância destes efeitos é função dos dois parã- 
metros adimensionais seguintes (ALMEIDA, 1981): 


t 
Cr = EA, Número de Courant ) (18) 
Á x 
e 
NO = E ( Número de onda ) (19) 
3 


O número de Courant permite uma quantificação do 
afastamento dos pontos de cálculo E e D relativamente 
aos nós À e B da malha de cálculo (Figura 3) e o número 
de onda dá uma medida do grau de discretização relativa- 
mente ao comprimento de onda característico. 

No caso do Cr ser igual à unidade não há interpolações 
e, por consequência, não ocorrem as alterações referidas. 
No caso de Cr ser inferior à unidade os efeitos em causa 
serão tanto mais significativos quanto maior for o valor de 
NO. Como regra geral, não são de aconselhar valores 
de Cr inferiores a 0,5 e valores de NO superiores a 0,1. Em 
valores absolutos, o amortecimento numérico vai também 
depender das caracteristicas do sistema elevatório e das 
características das perturbações que ocorrem durante o re- 
gime transitório. À caracterização destes dois factores físi- 
cos será feita, na parte II do presente artigo. por parâme- 
tros adimensionais. 


6. SAÍDA DE RESULTADOS 


O programa RESBOM fornece como saida (“output”) os 
valores das seguintes grandezas: 


— cota piezométrica, altura piezométrica e caudal, em 
cada secção de cálculo; 

— velocidade de rotação do grupo electrobomba; 

— envolventes das cotas piezométricas e das alturas pie- 
zométricas externas ao longo da conduta; 

— instante de ocorrência do fechamento das válvulas de 
retenção; 
caudal nas bombas no último instante de cálculo ou 
de simulação (controlo de que o tempo total de simu- 
lação é suficiente); 

— Intervalo de tempo teórico de funcionamento de um 
eventual “by-pass” no grupo electrobomba; 

— dados iniciais e parâmetros característicos (e.g. pe A); 


Quando o programa é utilizado num estudo sistemático, 
com variação dos valores dos parâmetros p e À, é também 
impresso um quadro com as variações extremas relativas 
das cotas piezométricas com todas as secções de cálculo. 


7. EXEMPLO DE UTILIZAÇÃO 


Uma das razões que poderá justificar a colocação dos 
grupos electrobomba numa secção intermédia da conduta 
elevatória é tornar mais económica a protecção contra o 
golpe de aríete ou mesmo evitar a colocação de dispositivos 
especiais de protecção. 

Na fase de análise preliminar, pode mostrar-se vantajoso 
estudar os efeitos, nas cotas piezométricas extremas, de 
diferentes localizações dos grupos electrobomba. A título 
de exemplo, considere-se um sistema elevatório com as se- 
guintes características gerais (Figura 4): 


Comprimento total 
Diâmetro da conduta D = 10m 


Celeridade das ondas elásticas c =800m/s 
Caudal inicial em regime permanente Q,= Imº/s 
Altura de elevação topográfica Hg = 25,6 m 
Velocidade inicial dos grupos No = 1000 r.p.m. 
Rendimento dos grupos no = 0,80 


Coeficiente da curva característica do grupo electrobomba 


A = 0,5037 x 104 m (r.p.m.)? 
= —0,1852 x 102 m (rpm)! (m/s)! 
C=0,1852x 102 m (m?/s)-? 


Perda de carga total inicial AH,=4m 
Parâmetro característico da conduta pn = 3,4 
Coeficiente de compensação energética A = 0,2 


CONFIGURAÇÃO & 


A 
A 
ao ga 


CONFIGURAÇÃO E 


FIGURA 4 — EXEMPLO DE APLICAÇÃO CONFIGURAÇÕES ESTUDADAS 


Consideram-se para comparação três configurações (Fi- 
gura 4) consoante o perfil da conduta e a posição do grupo 
electrobomba equivalente: configuração A com L|=0 e 
configuração Be € com L, = 750 m. As simulações foram 
efectuadas com Ax = 2580m e At = 03]s (Cr=1e 
NO = 0,10). No Quadro 7.1 indicam-se, para cada confi- 
guração e para cada secção de cálculo, os valores da cota 
do eixo (z) e das alturas piezométricas minima (y,) e má- 
xima (ym). 

Os resultados obtidos permitem concluir que a configu- 
ração A corresponde à ocorrência de pressões subatmosfé- 
ricas ao longo de toda a conduta com especial predominân- 
cia nas secções mais afastadas do grupo electrobomba, 
Mantendo o perfil da conduta (perfil linearmente ascen- 
dente), a colocação do grupo a 750 m do reservatório (con- 
figuração B) agrava substancialmente os valores das pres- 
sões mínimas na conduta de aspiração havendo, contudo, 
um benefício nas secções para jusante do grupo electro- 
bomba, exceptuando-se a secção S5 (distante 1500 m do 
reservatório de jusante) na qual há um agravamento. No 
que respeita as pressões máximas, a configuração B é mais 
favorável em todas as secções de cálculo. A configura- 
ção €, não obstante conduzir a um agravamento das pres- 
sões máximas na conduta a jusante do grupo electrobomba 
conduz a alturas piezométricas minimas sempre positivas. 

Deste exemplo pode concluir-se ser dificil resolver o pro- 
blema da ocorrência de pressões mínimas muito baixas 
com a deslocação do grupo electrobomba. Esta deslocação 
pode ser mais proveitosa no caso de o perfil da conduta ser 
favorável (perfil côncavo) e se o problema a resolver for o 
de diminuir as pressões máximas: na secção S5, por exem- 
plo, a altura piezométrica máxima para o perfil da configu- 
ração A é 36,6 m enquanto que para o perfil da configura- 
ção Bé 26,7 m. 

Para a configuração A seria mais eficaz aumentar o valor 
do PD” por intermédio de volantes de inércia pois deixa- 
riam de ocorrer pressões subatmosféricas e diminuiriam 
ligeiramente as pressões máximas, 

É de referir que na escolha do comprimento da conduta 
de aspiração há que atender às condições de funcionamento 
em regime permanente, em particular no que concerne a 
ocorrência da cavitação na conduta e nas bombas. 
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QUADRO 7.1 - Exemplo de aplicação. Cotas topograficas e valores das alturas 
piezomêtricas extremas para as tres configurações consideradas 
na Figura 4 (p=3,4 eàr=õo0,21) 
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PARTE II — Exploração do Modelo RESBOM Altura total de elevação inicial H ho 
Perda carga total inicial SH, 

8. PARÂMETROS ADIMENSIONAIS BÁSICOS Celeridade das ondas elásticas C 
Massa volúmica de água p 


Na exploração sistemática de um modelo, matemático 
ou físico, há que ter em conta a análise dimensional, tendo 
em vista a caracterização dos sistemas e das condições de 
cálculo por parâmetros adimensionais. A utilização ade- 
quada destes parâmetros permite uma generalização racio- 
nal dos resultados obtidos por via das simulações numéri- 
cas e uma apreciável redução no número de cálculos a efec- 
tuar. 

No estabelecimento dos parâmetros adimensionais há 
que distinguir as grandezas do fenómeno hidráulico (golpe 
de ariete provocado pela saída de serviço dos grupos elec- 
trobomba) e as grandezas que caracterizam os efeitos nu- 
méricos relacionados com a técnica de integração adop- 
tada. No que respeita o caso em estudo, as grandezas do 
primeiro grupo mais relevantes de acordo com as hipóteses 
básicas admitidas. são as seguintes (Figura |): 


Aceleração da gravidade g 


De entre as 12 grandezas intervenientes, escolhem-se três 
grandezas-base dimensionalmente independentes: p, Va € 
Hyo- De acordo com o teorema de Vaschy-Buckingham, o 
fenômeno hidráulico pode então ser caracterizado por nove 
parâmetros adimensionais (Ay, Ac. Ap, Ac Ago Ago Ag À 
e Ay) definidos no Anexo 2. À experiência na análise do 
golpe de ariete indica, contudo, que as condições de seme- 
lhança podem ser menos exigentes, sendo admissível consi- 
derar um menor número de parâmetros adimensionais por 
multiplicação de alguns dos parâmetros indicados no 
Anexo 2. Assim, os parâmetros adimensionais fundamen- 
tais escolhidos para o presente estudo são os seguintes 
(ALMEIDA e CUNHA, 1983): 


Variação de cota piezométrica AH 

Velocidade inicial do escoamento Vo Parâmetro característico da conduta, 
Velocidade de rotação inicial do grupo Ve 

Comprimento da conduta de aspiração L;, N ev 

Comprimento total da conduta elevatória Ls p= E A i 

Área da secção transversal da conduta elevatória Ss ; a 

Inércia das massas girantes (MDº - PD £) MD? 
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Coeficiente de compensação energética, 


iá PD? Nº 
A = comam neo 
N 2 

Y SL, Vê 


2 
A, Ay 
AA, 


(sendo V, = Ky No, com V, emrad/s e N9 em r.p.m.) 
L, L; 


AB GE” E ÀA3= 


Parâmetros de comprimento, 


Parâmetros da variação de cota piezométrica, Ay = 
(CAmm E Amgm) 


AH, 
H bo 


Parâmetro da perda de carga, À À = 


sendo Kw um parâmetro de conversão numérica de unida- 
dese A H as variações extremas de cota piezométrica (dife- 
rença entre as cotas piezométricas máxima e mínima e a cota 
piezométrica estática, respectivamente para A gym € A gm). 

O parâmetro característico p traduz a relação entre a 
energia cinética inicial da coluna liquida e a energia poten- 
cial de deformabilidade elástica. O coeficiente de compen- 
sação energética À traduz a relação entre a energia inicial 
das massas girantes do grupo electrobomba equivalente e a 
energia cinética inicial da coluna líquida. Desta forma, À 
pode ser definido do seguinte modo: 


p= — 20) 
E, (20) 
sendo: 
PD? Nº 
E, =—— (21) 
7153 
E; = 0,5 p SL; Vê (22) 


com PD em N.m?, No em r.p.m., p em kg/mº, S em 
m?, Ly em m e V; em m/s. 

No que respeita os parâmetros referentes aos efeitos nu- 
méricos foram seleccionados, para o presente estudo, os 
parâmetros C, (número de Courant) e NO (número de 
onda) que foram definidos em 5 (Parte 1). 

No que respeita as características hidromecânicas das 
bombas hidráulicas, a respectiva caracterização, tendo em 
vista a análise do golpe de ariete, pode ser tentada a partir 
do parâmetro N,, número especifico de rotações. 


9. METODOLOGIA SEGUIDA NA EXPLORAÇÃO 
DO MODELO 


Na exploração sistemática do modelo matemático RES- 
BOM (descrito na Parte 1) fizeram-se variar os valores dos 
parâmetros Ag, À € p: 


— consideraram-se oito posições para o grupo electro- 
bomba equivalente (valores de Ap compreendidos 
entre O e 0,7); 

— para cada valor de À p consideraram-se 11 valores para 
o PD? sendo o valor mínimo 300 N.m? e os restantes 
obtidos por sucessivas duplicações (valores de À 
compreendidos entre 0,025 e 26,5); 

— para cada par de valores de Ag e A fez-se variar o 
valor da celeridade entre 300 m/s e 1000 m/s o que 
corresponde a considerar oito valores para o parâmetro 
p (valores de p compreendidos entre 1,3 e 4,3). 


As simulações foram efectuadas com As = NO = 0,10, 
variando os valores de €, entre a unidade e 0,50 à excepção 
de duas situações correspondentes a À = 0,025 sendo, res- 
pectivamente, p = 1,3 (€C,= 0,34) e p= 1,7 (C, = 0,45). 
As variações dos parâmetros são efectuadas de modo 
automático de acordo com as instruções de um programa 
adicional do programa RESBOM. Foram também execu- 
tados alguns cálculos de sensibilidade para avaliação dos 
efeitos, nos resultados, de diferentes valores do passo de 
cálculo At (influência de C,) e da perda de carga inicial 
(influência de À À). No âmbito do presente trabalho execu- 
taram-se cerca de 1500 simulações numéricas tendo sido 
escolhido, como sistema elevatório base, o sistema com a 
configuração B do exemplo de aplicação apresentado em 7 
(parte 1) e considerando um grupo electrobomba equiva- 
lente com N, = 78 r.p.m.. 


10. RESULTADOS DAS SIMULAÇÕES 
10.1 — Objectivos do estudo 


No estudo efectuado, a exploração do programa de cál- 
culo ou modelo matemático RESBOM teve por objectivo 
principal a obtenção das variações extremas da cota pie- 
zométrica ao longo da conduta elevatória, para diferentes 
valores de Ag, À e p, decorrentes da saida de serviço do 
grupo electrobomba. Na continuação de estudos anteriores 
(ALMEIDA, 1981), houve também o interesse em obter o 
tempo de anulação do caudal na bomba (T,) e a velocidade 
de rotação do grupo no instante de fechamento da válvula 
de retenção (N,). No presente artigo apresenta-se uma sú- 
mula dos resultados obtidos e da análise efectuada os quais 
são descritos com pormenor em ALMEIDA e CUNHA, 
1983. 


10.2 — Variações das cotas piezométricas 


Às variações extremas da cota piezométrica são apresen- 
tados de forma adimensional pelos parâmetros À y cujos 
valores, em cada secção de cálculo, são obtidos do seguinte 
modo (Figura 1): 


— Secções a jusante do grupo electrobomba 


Aum = AHy/ Ho = (Hy — Zy)/ Ho (23) 
À Hm = À Hm [ Hho = (Hm — Z3) [ Hbo 
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— Secções a montante do grupo electrobomba 


Aum = à Hy/ Hi = (Hu — Zy)/ Hpo 


| (24) 
A Hm = A Hm / Hbo E (Hm e Zm) / Hho 


em que Hy = cota piezométrica máxima, H,, = cota 
piezométrica minima, Zm = cota do nível de água no re- 
servatório de montante, Z; = cota de nivel de água no 
reservatório de jusante e Hp, = altura total de elevação 
inicial, Às simulações foram efectuadas com 12 secções de 
cálculo correspondendo as secções 0 e 1 | respectivamente à 
secção de entrada na conduta de aspiração (junto ao reser- 
vatório de montante) e à secção de saida da conduta eleva- 
tória (junto ao reservatório de jusante). Nestas duas secções 
os valores de À y são sempre nulos. 

A exploração sistemática do programa de cálculo RES- 
BOM correspondeu à obtenção de um elevado número de 
resultados cuja compilação, apresentação e análise não é 
fácil de apresentar. No presente artigo faz-se unicamente 
uma sintese dos resultados mais importantes que poderão, 
eventualmente, ter interesse prático no projecto de sistemas 
elevatórios. 

Nas Figuras 5 e 6 indicam-se os valores de A ym e de 
À Hm: respectivamente a montante e a jusante do grupo, para 
cinco valores de À (A = 0,4; 0,8; 1,6: 3,3 e 6,6) e para di- 
ferentes valores de A p (diferentes localizações do grupo 
electrobomba equivalente) mantendo-se constante o valor de 
p (o = 3,0). Por seu turno, nas Figuras 8 e 9 indicam-se os 
valores de A ym e de A gm nas mesmas secções de cálculo 
anteriores mas em função do parâmetro p, para três valores 
de A (A = 0,4; 1,6€ 6,6) e para três valores de À p (A p= 0,2; 
0.3 e 0,4). Na Figura 7 indicam-se os valores extremos de 
Aum e de A gm do longo da conduta (envolventes) no caso 
de Ag = 0,3, paratrês valoresde À (A = 0,4; I6 e 6.6)e 
dois valores de p (p = 13 e 43). 

Dos resultados obtidos no estudo podem retirar-se, entre 
outras, as seguintes conclusões: 


para valores de À superiores a três, os efeitos elásticos 
deixam de ser significativos (influência pouco signifi- 
cativa do parâmetro p) passando as envolventes dos 
valores extremos de Ay a serem lineares (compor- 
tamento como coluna rigida); 


para valores de À inferiores a três, a influência de p 
[faz-se sentir de forma oscilatória (Figura 8 e 9) nos va- 
lores de À gym a jusante do grupo e de À gm à mon- 
tante do grupo; 

— em termos práticos e para valores de p superiores a 
dois, os efeitos elásticos deixam de ter influência rele- 
vanteem Ay para À > 0,8€cas envolventes passam 
a ser quase lineares para À > Lé6: 

— para elevados valores de À, os valores de À Hm cor- 
respondentes às mesmas secções de cálculo a jusante 
do grupo tendem a ser os mesmos qualquer que seja o 
valor de À p; 

— a ocorrência de cotas piezométricas no grupo possibi- 
litando o funcionamento de um “by-pass”, para pro- 
tecção contra o golpe de ariete, é muito difícil (e.g. 
paraA g = Otal só ocorre para p < 0,21). 
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Do ponto de vista prático, foi possível encontrar uma 
regra de redistribuição de valores de Ay nas secções 
adjacentes ao grupo, válida para elevados valores de 
A (A > 2,5Jo valor de A ym na secção de jusante do 
grupo no caso de Ap = O é igual à soma, em valor abso- 
luto, dos valores de À gym € de À Hm correspondentes, res- 
pectivamente, às secções de jusante e de montante do grupo 
no caso de Ag É O e, de modo equivalente, o valor 
de A Hm na secção de jusante do grupo para Ag = 0 é 
igual à soma, em valor absoluto, dos valores de A ym € 
de A Wm correspondentes, respectivamente, às secções de 
montante e de jusante do grupo no caso de Ag * 0. 


FIGURA 5- 
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FIGURA 6 = VALORES DE Jem E hmm NA SECCAO IMEDIATAMENTE A 
JUSANTE DO GRUPO PARA DIFERENTES VALORES DE 


10.3 — Tempo de anulação do caudal 


No Quadro 10.1 apresentam-se os valores de razão T,/ Ty. 
tempo de anulação do caudal no grupo electrobomba sobre 
tempo de inércia rigida, em função de alguns valores de À e 
de Ap. para p = 3,0. 


De acordo com os resultados obtidos, T, tende a ser 
superior a T, e pode depender de p. A influência deste 
último parâmetro tende a atenuar-se à medida que o valor 
de À aumenta: para Ag = 0, Ty / Ty varia entre 1,71 e 
1,49 e entre 5,06 e 4,60, respectivamente para valores de À 
de 0,4 e 6,6, correspondendo os limites de variação a valo- 
resde p de l3ec4,3. 


10.4 — Energia útil de rotação das massas girantes 


Contorme é conhecido (ALMEIDA, 1981), a energia ci- 
nética inicial das massas girantes de um grupo electro- 
bomba não é, em geral, totalmente cedida à coluna liquida 
após a anulação do binário motor e início do regime transi- 
tório. Com efeito, para valores elevados de A, o valor da 
velocidade de rotação do grupo no instante em que se anula 
o caudal na bomba, N,, tende a ser uma fracção significa- 
tiva do valor inicial N,,. 

No quadro 10.2 apresentam-se os valores da razão N,/ N, 

em função de alguns valores de À e Ap. para p= 1,3 e 
p = 4,3. Os resultados apresentados evidenciam que o 
comportamento como coluna rigida, isto é, 
a não dependência relativamente aos parâmetros p e Ap 
mantendo-se A constante, é admissível para valores de 
À > 33,nocasode p= Ijede A > 0,4nocaso de p = 4,3. 
Pode assim admitir-se que a energia útil das massas giran- 
tes E,, é a seguinte: 


Es =(1-Na/ No) E, (25) 


O valor de E, e, por consequência, a eficiência de um 
eventual volante de inércia ou de massas adicionais, para 
protecção contra o golpe de aríete, vai diminuindo à me- 
dida que o valor de À aumenta. 


ha E h(p=30 E ».=0,10) 


11. BREVE ANÁLISE DE SENSIBILIDADE 


Para avaliação dos efeitos nos resultados dos valores de 
At e da perda de carga total inicial AH, executaram-se 
alguns cálculos de sensibilidade com novos valores de €, 
(um terço dos anteriores) e com novos valores de 
A sn thÃ À = 0€ 0,14). Dos resultados obtidos podem reti- 
rar-se as seguintes conclusões: 


— nas secções imediatamente a jusante do grupo elec- 
trobomba, o valor de À ym é Um pouco afectado pela 
diminuição de At; pelo contrário, nas secções mais 
próximas do reservatório de jusante, os valores de À um 
mostram-se ser muito sensíveis à variação de At; 
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QUADRO 10.1 — Valores de 1, / Ty em função de A ge deA no caso dep =3,0(1,-0,5p Tp). 
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QUADRO 10.2 — Valores de N, / N; em função de À e de 
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— para valores crescentes de p acentua-se a diminuição 
dos valores de A yum motivada por amortecimento 
numérico (valores de €, menores); 

— a perda de carga faz diminuir, em valor absoluto, o 
valor de À Hm € diminuir as cotas piezométricas má- 
ximas; para elevados valores de A, as cotas prezomé- 
tricas não excedem as cotas correspondentes ao re- 
gime permanente inicial (À x = 0,14); 

— a perda de carga tende a retardar a anulação do cau- 
dal nas bombas e a diminuir o valor de N4 / No. 


Concluiu-se também que a determinação das cotas pie- 
zométricas mínimas considerando nulas as perdas de carga 
e admitindo que a cota do reservatório de jusante se eleva 
de AH,, o que equivaleria a manter o valor de Hy,, pode 
conduzir a resultados contra a segurança. 


12. CONSIDERAÇÕES E RECOMENDAÇÕES FINAIS 


O programa de cálculo RESBOM, desenvolvido no CE- 
HIDRO no âmbito do presente trabalho, está apto a execu- 
tar simulações sistemáticas tendo em vista a elaboração de 
tabelas de cálculo para análise dos efeitos do golpe de 
ariete. Os parâmetros adimensionais seleccionados e utili- 
zados no presente trabalho parecem ser adequados à carac- 
terização dos regimes transitórios em condutas elevatórias. 

Há ainda que esclarecer ou analisar com maior cuidado 
os efeitos numéricos (atenuação ou dispersão) dependentes, 


fundamentalmente dos parâmetros adimensionais designa- 
dos por números de Courant e de onda. É ainda de aconse- 
lhar a execução de uma análise sistemática dos efeitos re- 
sultantes das perdas de carga, da variação do rendimento 
dos grupos e de diferentes valores do parâmetro N, das 
bombas hidráulicas, A interacção das válvulas de retenção 
como escoamento e os eventuais atrasos ou antecipações no 
fechamento das mesmas poderão modificar os resultados 
obtidos e apresentados no presente artigo. À ocorrência de 
pressões subatmosféricas pode modificar o valor da celeri- 
dade e provocar a rotura da coluna líquida em zonas locali- 
zadas da conduta elevatória. Estes efeitos resultantes do 
fenómeno da macrocavitação transitória deverão ser consi- 
derados na análise dos regimes transitórios em sistemas 
hidráulicos reais (ALMEIDA, 1981). 
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ANEXO 1 
Os parâmetros adimensionais p e A, utilizados neste 
artigo e que são apresentados na Parte II, são definidos do 
seguinte modo: 
à-= Er | Ec 


sendo 


Eg = PD” Nº | 7153 (energia cinética das massas girantes) 


Ec = D35p E; 8 q (energia cinética da coluna líquida) 


ae] 


p = 6 Vo | (£Hbo) 


em que PD» = parâmetro característico da inércia do grupo 
electrobomba em N.m?; N, = velocidade de rotação inicial 
em r.p.m.; Ly = comprimento total da conduta; S = área da 
secção transversal da conduta; V, = velocidade inicial do 
escoamento em m/s; c = celeridade das ondas elásticas em 
m/se Hy, = altura total de elevação das bombas em mc.a.. 
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ANEXO 2 


Os nove parâmetros adimensionais referidos no presente 
artigo são os seguintes: 


|) Parâmetro da inércia do grupo 


MD? PD? 
É “o 5 
p Ho É 4 Ho 
S 
21 Parâmetro da secção da conduta As= 
Ho 
PAN q e 
3) Parâmetro da perda de carga Mi a 
bo 
c 
4) Parâmetro da celeridade (inverso do Aç-= 
número de Mach) Vo 
| £ Hy 
5) Parâmetro da gravidade Ag — - 
“a 
L, 
6) Parâmetros de comprimentos AB= E 
3 
L; 
Aa = | 
Hho 
Ve Hho 
7) Parâmetro da velocidade periférica  Ay= 
do grupo Vo 
AH 
8) Parâmetro de variação de cota Ay = 
piezométrica Hho 
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INTERATÔMICAS A PARTIR DA ESTRUTURA 
FINA QUE SE PROLONGA PARA ALÉM DA 
DESCONTINUIDADE DE ABSORÇÃO DOS RAIOS X 


MARIA ISABEL GASPAR DE BARROS ALVES MARQUES 
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RESUMO 


É projectada uma montagem experimental que permite 
obter um espectro de EXAFS. A montagem descrita utiliza 
como fonte o radiamento de distribuição espectral continua 
duma ampola de raios X, de ânodo rotativo. Os espectros 
obtidos, por transmissão, permitem a determinação de dis- 
tâncias interatômicas com um rigor de 0,01 A. É feita uma 
breve referência à análise de Fourier da função x (k), que 
dá conta das variações relativas que sofre o coeficiente de 
absorção ao longo do espectro: conhecidos os parâmetros 
intratômicos, as diferenças de fase e a amplitude da onda 
rectrodifundida, é possivel obter grandezas interatómicas 
— q distância entre dtomos vizinhos e o número de coorde- 
nação respectivo. 


1. INTRODUÇÃO 


A estrutura fina que se prolonga para além da desconti- 
nuidade do espectro de absorção dos raios X duma deter- 
minada substância, conhecida por EXAFS — “Extendend 
X-ray Absorption Fine Structure” — é devida às oscilações 
do coeficiente de absorção dos raios X que se podem obser- 
var quando os valores da energia variam desde 50 eV até 
cerca de 1000 eV acima da descontinuidade de absorção. 

A primeira teoria interpretativa do fenómeno, devida a 
Kronig, surgiu em 1931 para o caso dum cristal (estado 


Original recebido para publicação em 18/12/84. 


sólido) e um ano depois para um gás. O mau acordo quan- 
titativo entre as primeiras teorias que apareceram, € os re- 
sultados experimentais, fizeram com que só quarenta anos 
depois, em 1971, com a publicação das investigações de A. 
Stern, F. W, Lytle e E. Savers se pudesse utilizar a EXAFS 
como técnica experimental. 

A EXAFS permite a determinação de distâncias inter- 
atóúmicas na matéria condensada. A vantagem em relação à 
difracção dos raios X é permitir conhecer a estrutura em 
torno de cada espécie de átomo separadamente. 

Este método de investigação da estrutura interatómica 
tem sido na última década aplicado a sistemas de multi- 
componentes, na biologia, na química inorgânica, para de- 
terminar a estrutura de soluções aquosas, de catalisadores, 
ou detectar a presença de impurezas em sólidos. 

A EXAFS necessita duma fonte de raios X contínua. 
Esta técnica teve um grande desenvolvimento quando se 
pode usar como fonte a radiação de um sincrotrão, que 
permite obter os mesmos espectros durante menores tem- 
pos de exposição, 


2. O FENÔMENO DA EXAFS, SEU FORMALISMO 


O fenómeno da EXAFS, as oscilações do coeficiente de 
absorção dos raios X em função da energia, foi explicado 
por Kronig considerando a modificação do estado final do 
totoelectrão pelos átomos que circundam o átomo central, 
quando sobre este incide um fotão de raios X de energia 
superior à descontinuidade de absorção. A modificação do 
estado final foi atribuída ao efeito das interferências junto 
do núcleo resultantes da onda de matéria associada ao fo- 
toelectrão que parte do átomo absorvente e que é rectrodi- 
fundida por cada um dos átomos vizinhos. As interferên- 
cias dão origem aos máximos e mínimos no espectro da 
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EXAFS de acordo com a diferença de fase correspondente. 
As amplitudes das oscilações podem ser da ordem de 10% 
ou superiores em relação à curva monotónica a seguir à 
descontinuidade de absorção dum átomo isolado. 


A expressão que traduz este fenómeno pode escrever-se 


u—u N, 
jtig= =D e= E 6d e) 
Ho ] Fk 
2r 
e WS 
senQk nte 1 e Mb q) 


em que 

u é o coeficiente de absorção linear medido experimen- 
talmente 

ua é o coeficiente de absorção do átomo não perturbado 


é Ao A ? 

k =—— é a prandeza do vector número de ondas do 
À 

fotoelectrão 


|$; (k 7) | é à amplitude da onda rectrodifundida elasti- 
camente 


& (ngm 6 gy O (2) 


r; é a distância do átomo absorvente à camada de coorde- 
nação | 

N; é o número de átomos à distância r; 

2k rj no argumento do seno, traduz a diferença de fase 
introduzida pela ida e volta da onda desde o átomo central 
a um átomo vizinho 

YU (k)=60 (k) + 26, (k) éa diferença de fase total sofrida 
pelo fotoelectrão, soma da diferença de fase 6 (k), da onda 
rectrodifundida, com o dobro da diferença de fase ô, (k) 
introduzida pelo potencial do átomo absorvente (fotoexcitado) 


E e 
Iak* 
meo A é um factor do tipo de Debye-Waller que no ca- 


e 
so dum sólido faz intervir o efeito das vibrações térmicas € 
num liquido as variações de posição das partículas 

2 

» (K) é um factor de declínco que se torna necessário in- 
troduzir porque o efeito de interferência se atenua tanto mais 
quanto maior é a distância dos átomos vizinhos ao átomo 
central, é habitual chamar a À; (k) livre percurso médio 
electrão - electrão. 

As modulações em x (k) são assim determinadas por 
|f; (k, m)| e pelo termo sinusoidal, As outras contribuições 
são independentes de k ou decrescem monotonicamente 
com &. 


3. A MONTAGEM EXPERIMENTAL 


Para obter um espectro de EXAFS é necessário dispor 
duma fonte de raios X que forneça um espectro continuo 
intenso e dum espectrómetro com um adequado sistema de 
detecção. 
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A. A Fonte 


A fonte de raios X escolhida para este projecto é uma 
ampola de raios X de ânodo rotativo e de foco fino des- 
montável, que permite obter um espectro continuo suficien- 
temente intenso. 

O elemento do anticátodo pode ser a prata, o molibdênio 
ou outro de acordo com o espécimen em estudo, uma vez 
que os anticátodos são fáceis de substituir neste tipo de 
ampolas. 

O cátodo é constituído por um filamento helicoidal de 
tungsténio que permite fazer incidir sobre o ânodo o feixe 
de electrões numa região cujas dimensões se podem modifi- 
car por simples adaptação de uma cúpula metálica. Assim, 
para estudos da descontinuidade de absorção K de ele- 
mentos de número atómico médio, a “mancha” focal esco- 
lhida tem 0,5x5mm e para elementos de número atómico 
elevado 0,2x2mm. 

Este tipo de ampolas (de ânodo rotativo) pode apresen- 
tar por vezes problemas de instabilidade mecânica que se 
traduzem numa certa tremura da mancha focal sobre o 
ânodo. A ampola escolhida deverá ser de uma marca que 
garanta uma estabilidade mecânica de 5 um. 

O ângulo de saída do feixe pode variar entre O e 17. Es- 
colhemos um ângulo de saida da ordem de 3º o que dá para 
dimensões aparentes da mancha focal sobre o ânodo (pro- 
jectada segundo o ângulo de 3º) 0,025x5mm e 0,0]x2mm 
respectivamente. 

O ânodo tem a velocidade de rotação de 3000 ou 6000 
r.p.m. de acordo com a potência que pretendemos utilizar. 

A ampola pode trabalhar com uma tensão máxima de 
60kV e a corrente máxima de 300 mA, o que corresponde a 
uma potência máxima de 15 kW (para uma mancha focal 
de 0,5x10mm e para um anticátodo de cobre). Esta potên- 
cia é dez vezes superior à potência que pode suportar uma 
ampola de ânodo não rotativo para o mesmo anticátodo e 
para dimensões da região focal da mesma ordem de gran- 
deza. Para as dimensões das manchas focais escolhidas os 
valores da potência máxima que a ampola pode suportar 
são respectivamente de 10,0 kW e de 3,0 kW. 

A estabilização da tensão do cátodo e da corrente da 
ampola impede que as oscilações destas grandezas sejam 
superiores a 0,1%. 


Uma ampola deste tipo necessita, para refrigerar O 
ânodo, de 551 de água por minuto a uma temperatura infe- 
rora 15º Ce a uma pressão de 4 kg/cm2, o que exige um 
sistema de refrigeração em circuito fechado. 

Um sincrotrão produz um feixe com maior densidade de 
fotões (cerca de 10º vezes superior, tendo em conta já a 
grande distância da fonte ao espécimen) e isto é uma vanta- 
gem relativamente às ampolas de ânodo rotativo, mas estas 
dão feixes cujas características são mais estáveis no decor- 
rer do tempo, o que é muito importante. 

As ampolas fornecem além do espectro continuo o espec- 
tro característico do material que constitui o anticátodo (e 
o de uma ou outra impureza). Há pois que ter em conta esta 
circunstância escolhendo um anticátodo adequado ao espé- 
cimen a estudar, As riscas de emissão caracteristicas do 
anticátodo podem ser utilizadas para calibrar a escala das 
energias, como se verá adiante e em particular as riscas de 
emissão duma impureza como o tungsténio, por exemplo, 
podem servir para verificar a linearidade da escala das 
energias. 


A fonte escolhida permite obter o mesmo rigor que um 
sincrotrão na determinação de distâncias interatómicas, exi- 
gindo apenas maior tempo de exposição. Esta fonte tem 
ainda a vantagem de ser mais acessivel à maior parte dos 
laboratórios de investigação. 


B. O Espectrómetro 


O espectrómetro deste projecto será constituído por um 
goniómetro horizontal 8 - 8, no eixo do qual se monta um 
monocromador, e por um detector. Um goniómetro deste 
tipo permite que o monocromador rode de um ângulo 8 
enquanto o detector roda de um ângulo 298. 


a. O monocromador 


O monocromador é um cristal curvo montado de forma 
a reflectir o feixe proveniente da ampola. Outras opções 
seriam menos adequadas: um cristal plano reflecte com 
menor intensidade do que um cristal curvo; as montagens 
por transmissão permitem obter menor intensidade ainda 
que as montagens por reflexão. 
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Fig.l — Esquema da montagem experimental: C-cátodo, A-ânodo 
(rotativo), M-cristal monocromador com os planos recticulares 
curvados segundo 2R e a superficie exterior segundo R, F-fenda de 
bordos cilindricos, D-detector. O ponto O é o centro do circulo de 
Rowland de raio R,e MF coincide com o raio do goniómetro. a 
representa o ângulo de saída do feixe de raios X e 6 é o ângulo de 
Bragg. 


A prensa deverá ser magnética para garantir a necessária 
regularidade da curvatura do cristal, 

As dimensões úteis da lâmina cristalina (depois de mon- 
tada na prensa) são 55xêmm. 

A geometria escolhida foi a de Johanson (4). Assim as 
lâminas são talhadas com uma curvatura 2R e encurvadas 
numa prensa de raio R. 


b. O Goniómetro 


No eixo do goniómetro é montada a prensa com o cristal 
monocromador. O goniómetro que escolhemos com um raio 
de 20cm deverá assentar sobre uma plataforma que permita 
translacções paralelas à direcção do feixe de raios X que sai 
da ampola. A distância AM do centro do anticátodo ao 
monocromador será assim regulável, Por sua vez o sistema 
detector e fenda deverá também ser regulável, segundo a 
direcção do raio do goniómetro, para permitir ajustar a 
posição de focagem. Isto pode fazer-se substituindo o su- 
porte fixo habitual do detector dum goniómetro horizon- 
tal, por um suporte munido dum “granzep”, que permite 
que o detector execute translacções finas segundo o raio do 
goniômetro, 

Tendo sido escolhida para o monocromador a geometria 
de Johanson, o centro do anticátodo da ampola de raios X 
assim como a fenda que precede o detector devem ser colo- 
cados sobre a circunfereficia (de focagem) que limita o cir- 
culo de Rowland (2) de raio R, fig. 1. A distância AM é 
igual à distância MF (do centro do monocromador ao 
ponto de focagem sobre a fenda do detector). O cristal 
monocromador terá assim os planos recticulares usados 
paralelos à superficie de talhe da lâmina do cristal (o = 0). 


A natureza e características do cristal monocromador 
deverão ser escolhidas de acordo com a energia da descon- 
tinuidade de absorção do elemento em questão, Assim, 
para estudar o espectro de EXAFS que se estende a seguir à 
descontinuidade K de um elemento de número atómico 
médio, como por exemplo o bromo (Z = 35), a lâmina do 
cristal monocromador que se sugere é o fluoreto de lítio 
talhada paralelamente aos planos (200). Esta lâmina tem 
um poder reflector intenso, e a distância d entre os planos 
recticulares vale 2d = 4,028 A. O ângulo de Bragg para a 
descontinuidade de absorção K do bromo, à qual corres- 
ponde a energia de 13473,7 eV é 9 = 13,21”. Para uma 
prensa de raio R = 450mm obtém-se AM = MF = 200mm; 
este cálculo foi feito não para o valor de 8 indicado mas 
para um valor 6 4 = 12,85' compreendido entre o maior 
e o menor valor do ângulo de Bragg, usado para se descre- 
ver um espectro de EXAFS. O raio do circulo de Rowland 
foi escolhido de forma a que fosse possivel obter um espec- 
tro de EXAFS intenso com uma resolução suficiente. 

Para as distância referidas, e utilizando a mancha focal 
com a maior dimensão na posição vertical, a janela de en- 
trada do monocromador, que se encontra a 160mm do an- 
ticátodo, deverá ter uma largura de $&mm e uma altura de 
5.5mm. Estas condições produzirão sobre a lâmina do mo- 
nocromador uma mancha cujas dimensões máximas são 
4Smm de largura e 6mm de altura, o que é aconselhável, em 
face das dimensões úteis da lâmina (55x8mm). Evitar-se-á 
assim que o feixe incida nos bordos da prensa ou em re- 
giões limitrofes da lâmina cristalina. A janela de saída do 
monocromador, a 40mm do seu centro terá uma largura € 
uma altura ligeiramente superiores respectivamente a &mm 
e a 6mm, afim de evitar que o feixe que é reflectido pelo 
monocromador seja diafragmado pelos bordos da janela de 
saida; esta permite evitar a passagem de raios X parasitas, 
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Embora a mancha focal sobre o anticátodo seja usada 
com a maior dimensão vertical, a sua altura e a geometria 
da montagem não exigem a utilização de fendas de Soller. 
À dimensão do ânodo projectado de 0,025xSmm corres- 
ponde sobre a fenda de entrada do detector uma mancha de 
altura inferior a 7mm o que dá uma divergência em altura 
de 0,1”, 

A fenda do detector deverá ter a largura de 0,025mm. 

Entre a janela de saída da ampola e a janela de entrada 
do monocromador será montado um colimador cilindrico; 
e outro entre a janela de saída do monocromador e a fenda 
de entrada do detector. 

Para valores da energia desde 200 eV abaixo até 1000 eV 
acima da energia da descontinuidade K do bromo (13473,7 
cV) o detector terá que percorrer uma região 20 = 2,26 .0 
intervalo em comprimentos de onda para o qual teremos 
que fazer leituras variará assim entre Amin = 0,8566 A e 
Amax = 0,9340 À. O elemento mais adequado para anti- 
cátodo será a prata cujos comprimentos de onda das riscas 
características variam entre 0,5638 À e 0,4860 À. 

Para o estudo do espectro de EXAFS que se estende a 
seguir à descontinuidade K dum elemento de elevado nú- 
mero atómico, como por exemplo o índio (Z = 49),a mon- 
tagem terá que sofrer algumas alterações: o cristal mono- 
cromador, que deverá ter também um poder reflector in- 
tenso, terá uma constante de rede menor porque o valor do 
ângulo de Bragg é agora menor também. Assim escolhemos 
um cristal de fluoreto de lítio (220) para o qual 2d = 2,848 A. 
Afim de se usar o goniómetro com o mesmo raio, 20cm, 
como neste caso o ângulo de Bragg é igual a 8,96' (calculo 
para a descontinuidade de absorção K do In em que o valor 
da anergia é igual a 27939,9 eV) o cristal deverá ser talhado 
com um raio 2R = 1300mm e encurvado numa prensa de 
raio R = 650mm o que nos dará como anteriormente 
AM = MF = 200mm, 

As dimensões da janela de entrada do monocromador 
(usando também a maior dimensão da mancha focal na 
vertical) serão agora 6mm de largura por ômm de altura; a 
janela de saida deverá ter qualquer das duas dimensões 
ligeiramente superiores a 6mm. 

Como neste caso o ânodo projectado tem a dimensão de 
0,0]x2mm, a dimensão da mancha sobre a fenda do detec- 
tor é de 9mm e a sua divergência é de meio grau. Não se 
justifica também neste caso a montagem de uma fenda de 
Soller; se a dimensão e a posição da janela de entrada do 
detector o exigissem poderia usar-se uma fenda que limi- 
tasse em altura, 

A fenda do detector deverá neste caso ter a largura de 
0,0 mm. 

Para aumentar a intensidade de um espectro de EXAFS 
é possível introduzir ainda um melhoramento de acordo 
com (42) que consiste em tornar os colimadores e a caixa 
do monocromador estanques enchendo-os então de hélio. 
Evita-se assim a absorção pelo ar ao longo do trajecto do 
feixe. Uma hipótese menos dispendiosa seria fazer vácuo 
não permitindo que a pressão desça abaixo duns quinze a 
vinte milimetros de Hg. Esta limitação é sugerida por causa 
de se ter que manter rigorosamente a composição química 
da amostra durante a experiência, que se situa no trajecto 
do colimador entre o monocromador e a fenda de entrada 
do detector. 

O anticátodo da ampola deverá ser de molibdénio uma 
vez que as riscas K características do anticátodo de prata, 
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que têm comprimentos de onda compreendidos entre 
0,5638 À e 0,4860 À , ficariam muito próximos do com- 
primento de onda A max = 0,4453 À do espectro de EXAFS 
do índio. Uma estimativa feita a partir do cálculo da lar- 
gura À À do dublete K a da prata quando a tensão apli- 
cada à ampola é da ordem de 35KV deu 0,09 À para essa 
proximidade (os valores usados para este cálculo foram 
extraídos de (1) ). 


C. O Detector 


A medição do coeficiente de absorção em função do va- 
lor da energia para espécimens concentrados pode fazer-se 
directamente medindo a atenuação do feixe incidente numa 
experiência de transmissão. Para este fim será adequado 
usar uma câmara de ionização, a não ser que a relação sinal 
— ruido seja pequena; então convirá usar antes um detec- 
tor de cintilação, por exemplo para soluções aquosas con- 
centradas. 

Os detectores de cintilação são adequados para radiações 
de comprimento de onda compreendidos entre 0,2 e 2,5 A. 
Têm uma resposta bastante linear e o seu tempo morto é da 
ordem dos décimos de microsegundo. 

A eficiência de contagem destes detectores é elevada, par- 
ticularmente para comprimentos de onda menores que | 
(é o nosso caso) o que reduz o tempo de acumulação de um 
número préfixado de impulsos para cada medida. O valor 
da capacidade do detector é da ordem dos 200 000 imp. /s. 


D. O Espécimen 


Devem adoptar-se duas precauções importantes: conse- 
guir um espécimen homogéneo e manter rigorosamente a 
sua composição química durante a experiência. 

Se se trata dum sólido será reduzido a pó e pode introdu- 
zir-se entre duas fitas adesivas que o protegem da humi- 
dade. Se se trata dum liquido terá que ser introduzido 
numa célula dum material quimicamente inerte (em relação 
ao espécimen pelo menos) com duas janelas plano-paralelas 
duma substância pouco absorvente, por exemplo “kapton” 
ou colódio (6). 

O espécimen pode ser colocado antes ou depois do mo- 
nocromador. Se for colocado depois evita-se a possibili- 
dade de se focarem os radiamentos de fluorescência que 
então se sobreporiam ao espectro de EXAFS, 


E. Crítica 


A resolução em energia da montagem pode obter-se 
a partir da expressão que dá o valor da energia E, em 
função do comprimento de onda A, e da lei de Bragg 


E - he he j 
À 2dsenô (3) 


em que h designa, como habitualmente a constante de 
Planck e c a velocidade da luz 


Calculando A E: 


AE = E cotg6 AB (4) 


obtém-se, para valor mínimo do incremento À 8 da mon- 
tagem, o valor do incremento que sofre o valor da energia 
duma observação para a imediatamente seguinte. 

Para o espectro de EXAFS que se estende a seguir à 
descontinuidade de absorção K de um elemento de nú- 
mero atómico médio (Z = 35) obtém-se na montagem 
projectada 


NE =7,2 eV 


e para um elemento de número atómico elevado (Z = 49) 
obtém-se 


SE = 89 eV 


Qualquer dos valores obtidos para a resolução nas duas 
montagens que projectamos é suficiente se atendermos a que 
a extensão espectral dos menores detalhes dum espectro de 
EXAFS é da ordem de 30 eV (cálculo feito para A k = 0,8 
A! e para um valor médio de k, por exemplo para 
k=5 A. 4 

Se a resolução da montagem não fosse suficiente podia 
melhorar-se reduzindo a dimensão da fenda do detector e 
consequentemente a dimensão do foco da ampola (ou pelo 
menos o ângulo de saída desta), ou ainda aumentando a dis- 
tância AM = MF. Qualquer destas duas hipóteses signi- 
ficaria porém uma melhor resolução à custa de uma menor 
intensidade. 

Os elementos sugeridos para anticátodos nos dois projec- 
tos são diferentes. Caso não fosse fácil substituí-los poderia 
usar-se sempre o molibdénio. 

Um detector de cintilação associado ao seu sistema elec- 
trónico permite normalmente rejeitar as harmónicas de or- 
dem superior, mas se não dispuséssemos dum sistema de 
detecção com esta possibilidade o tipo de ampola proposta 
permitiria eliminá-las satisfatoriamente, dadas as suas ca- 
racterísticas eléctricas. 

Para o exemplo dado do bromo pode ver-se no quadro | 
os comprimentos de onda das harmónicas de 2.º e 3º or- 
dem. Se admitirmos que esses comprimentos de onda calcu- 
lados são os valores a que corresponde a intensidade má- 
xima do espectro continuo, da expressão de (14) 


3 
— Amin (5) 


Max 2 


podemos calcular os valores de Amin € finalmente a lei de 
Duane e Hunt 


he 
eV = 


6 

À min (6) 

permite-nos calcular os correspondentes valores da tensão 
a aplicar à ampola. 


QUADRO 1 


Confronto entre o valor do comprimento de onda À (n) 
das harmónicas de ordem n, e os valores de Ain do 
espectro contínuo para diversos valores da tensão V. 


n À (n) e V 
(A) (À) (kV) 
| 0.92 0.61 20 
2 0.46 031 40 
E 031 0.20 62 


A tensão deverá ser igual ou superior a 20 kV. Se quiser- 
mos ter a certeza que não é produzido radiamento corres- 
pondente à harmónica de segunda ordem o valor da tensão 
a aplicar à ampola deverá ser inferior a 26,95 kV a que 
corresponde um Amin = 0,46 À, que é o comprimento de 
onda da harmónica de 2.º ordem. Como para o foco esco- 
lhido a potência máxima é de 10 KW, a corrente máxima 
que poderemos aplicar à ampola é de 300 mA para a tensão 
de 25 kV por exemplo. 


4. A OBTENÇÃO DUM ESPECTRO DE EXAFS 
A. A Medição das Energias 


Numa experiência de EXAFS pretende medir-se o coe- 
ficiente de absorção dos raios X em função da energia do 
feixe incidente. Assim obtém-se um espectro de emissão, 
fig. 2a, para energias que vão desde cerca de -200 eV a 
+ 1000 eV em relação à descontinuidade de absorção do 
elemento da substância para a qual se pretende obter o 
espectro de EXAFS. Com o espécimen no trajecto do feixe 
obtém-se na mesma região um “espectro de transmissão”, 
fig. 2b. 


Fig.2 — Valores da intensidade, em imp. /s, em função dos valores 
da energia, em eV, dos espectros que se obtêm quando o especimen 
está intercalado (1,), e quando não está (I,), no trajecto do feixe. 
Na região das baixas energias do espectro vê-se uma risca caracte- 
ristica e na região das altas energias uma risca duma impureza do 
anticátodo. 


O espectro de absorção obtém-se calculando pela lei de 
Lambert-Beer 


(7) 


em que ! é a espessura do espécimen e |, e 1, são respec- 
tivamente a intensidade que se obtém quando o espécimen 
está ou não intercalado no trajecto do feixe, fig. 3a. 


— 8 = 


Fig.3 — Espectro de absorção u (E) 1 
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As experiências que nos permitem obter os espectros de 
emissão e de transmissão dão-nos os valores das intensida- 
des para os valores duma certa gama do ângulo 2 8 lidos 
no espectrómetro. para converter cada valor de 9; em E, 
utiliza-se a expressão de (11) 


12398,52 
j 2 dsen 9, 


em que os valores da energia E, são dados em ev. 


Porém, para se obter os valores de E, com a precisão 
necessária (distâncias interatômicas r com a precisão de 
AO! À) a constante 2d a usar deverá ser, não a do forne- 
cedor da lâmina cristalina, mas a que se obtém aferindo-a 
por meio duma risca característica do espectro. Este pro- 
cesso tem a virtude de evitar as correcções de temperatura 
ou de polarização de Lorentz. 

Para os elementos bromo e índio os valores de 6. 


inic € O fin 
são os que figuram no quadro 2. 


QUADRO 2 
inte 9 in 
Bromo 13,41" 12,28" 
Índio 9,03 8,71º 


Por estes valores se vê que as descontinuidades K do 
bromo e do índio distam suficientemente entre si para os 
dois espectros de EXAFS se poderem obter sem sobreposi- 
ção. Poder-se-á assim investigar separadamente os dois espectros 
de EXAFS para estes dois elementos numa solução aquosa 
de In Bra. 

Para uma banda de energia de 7,2 eV um incremento da 
ordem de 3,6 eV é adequado. Isso corresponde, no caso do 
bromo,a À 26 = 0,014 ea 333 leituras e no caso do índio 
ai 26 = 0,006" e a 300 leituras, uma vez que A E = 89evV. 
Estes incrementos não levantam dificuldades experimen- 
tais. Dum modo geral um goniómetro permite incrementos 
mínimos de A 26 = 0,001", 


B. A Conversão da Escala das Energias 
numa Escola de Números de Onda 


Para se poder extrair dos resultados experimentais da 
EXAFS as informações relativas à estrutura da matéria que 
se investiga é necessário conhecer o coeficiente de absor- 
são u em função da grandeza k do vector número de on- 
das do fotoelectrão que é expulso do átomo perturbado. A 
função u (k) obtém-se de y (E) relacionando a expressão da 
energia cinética do fotoelectrão com a energética h v do 
fotão de raios X incidente 


h?k? | 
— = hv- E, (K)=E-E, (k) (9) 
im 


am |E-Eç(k)| | 
a (10) 
em que £, (k) é a energia de referência da descontinui- 
dade de absorção. 

A aplicação desta expressão tem a dificuldade de obter 
um valor de k com rigor, uma vez que o valor de E, é 
desconhecido. Rodeia-se esta dificuldade considerandoo 
um parâmetro de valor ajustável. Usa-se em geral um valor 
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de E, para o qual a curva que representa a descontinuidade 
de absorção tem uma certa particularidade tal como um 
ponto de inflexão. Existem na literatura valores de E, cal- 
culados teoricamente para o Nie o Nb (13). 


C. A Determinação de x (k) 


Para se obter x (k), partindo do coeficiente de absorção 
total qu (E), pode começar-se por subtrair a q (E)t a função 
que se obtém extrapolando os valores do coeficiente de 
absorção correspondentes à região lentamente variável que 
fica à esquerda da descontinuidade de absorção, fig. 3c à 
qual ajustamos os coeficientes dum polinómio de Victoreen 
(11) CAº + DA?. Obtém-se assim uma função uy (E) t, 
fig. 3b, devida somente à camada K do elemento em ques- 
tão. 

Para determinar x (k) precisamos ainda de conhecer o 
valor de u, (E), o coeficiente de absorção do átomo iso- 
lado, para o qual os cálculos teóricos não dão a precisão 
necessária, para os resultados que pretendemos obter 
(0,01%). Admitindo que os valores de uy (E) t coincidem 
com os que podem ser fornecidos a partir do declive suave 
da curva que representa up (k) t, a parte oscilatória res- 
tante é considerada À yu = ui (E) — ua (E) e obtém-se assim 


Ô u 
Ko 


x (k) = (11) 


normado para uma escala referida a um átomo, Para esta 
segunda fase ajusta-sea uy (k)t um polimónio em k, fig. 4. 


Fig.4 — Espectro de absorção u, (k)1 


Este processo que se encontra bastante usado na literatura 
pode de acordo com (41) conter erros sistemáticos que se 
poderão evitar usando polinómios “Splines” ou “B Splines” 
(12 e 26). 


Na fig. 5 pode ver-se representada uma função x (k). 


Fig.5 — Função x (k) obtida experimentalmente. 


5. DETERMINAÇÃO DAS DISTÂNCIAS 
INTERATÓMICAS 


A transformação de Fourier da função x (k) permite 
obter as distâncias interatómicas. 

A distância do átomo central aos átomos vizinhos mais 
próximos pode no entanto determinar-se por um método 
gráfico simples (23). Este método admite que a difusão de- 
vida à primeira camada predomina na função x (k). Assim, 
considerando apenas a ondulação principal, o periodo da 
função seno corresponderá ao da primeira camada de co- 
ordenação. A análise depende pois da variação da função 


sen [2k r; + 2 We; (k)] 
Se y; (k) varia linearmente com k 
W; (k) = -qjk+ B, (12) 


Substituindo no argumento da função seno, obtém-se 
para a primeira camada 


2k (ri; —-a|) +28/=5 17 (13) 


em que n = 0,2,4,... para os máximos da função x (k) e 
n = 1,3,5,... para os minimos. Um gráfico de n em função 
dos valores de k para os quais a função x (k) é máxima 
ou mínima permite obter r; — a, (do coeficiente angular 
da recta) e 8, (da ordenada para k = 0) 


Fig.6 — Valores de n, da expressão (13), para os diferentes valo- 
res de k do argumento do seno (expressão (11)) em que é máxima 
ou minima a função x (k). 


Utilizando este método para uma substância em que rj é 
conhecido determina-se o valor de a ,. Para uma substância 
de composição quimica semelhante, em que rj é desconhe- 
cido, usando o mesmo valor de a | pode obter-se a distância 
interatómica. 

A análise de Fourier da função x (k) pode fazer-se, em 
primeiro lugar, por meio do módulo da transformada com- 
plexa desta função. 

A transformada de Fourier de k”" x (k) no espaço dosr 
é dada por 


k 
E (1) = gts | DO qn x(k) W(k) e kr dk (14) 
k 


mun 


O factor k” destina-se a compensar o facto de a ampli- 
tude do espectro de EXAFS se reduzir quando o valor de 
cresce. Para uma região de k entre 4e 15 A"! os valores 


de n recomendados em função do número atômico dos 
atómos rectrodifusores figuram no quadro 3. 


QUADRO 3 


Valores de n usadas para expoente de k, na inversão de 
Fourier de x (k), relativos a diversos elementos, de 
número atômico Z. 


Elementos 


à ido Valores de n 
de número atômico Z 


L< M% É, 
6<Z< 5 2 
SI< Z< 86 | 


Nas fig. 7 e 8 pode observar-se o efeito de multiplicar 
x (k) por k” para o caso de um espectro de EXAFS do 
germânio (Z = 32). As amplitudes que se observam na fig. 8 
são obviamente mais uniformes que as que se observam na 
fig. 7, antes de x (k) ter sido multiplicado por k*. O efeito 
de multiplicar por k* começa a fazer-se sentir a partir de 
k =4 A"! que corresponde à energia dum fotoelectrão de 
cerca de 60 eV acima da descontinuidade de absorção. 


Fig.7 — Andamento da função x (k) tal como é obtida da expe- 
riência, 


Fig.8 — Andamento da função kº x (k) em que x (k) está repre- 
sentada na fig. 7, 
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Na expressão (14), W (k) é uma função com o perfil de 
Gauss, que tem o valor 0,1] em Kmin € EM Kmax (0). 

O módulo da transformada de Fourier complexa da fun- 
ção k”" y (k) W (k) a que chamaremos “função de estru- 
tura” é uma função radial do tipo da que representamos na 
hg. 9, Desta função obtêm-se as distâncias do átomo cen- 
tral aos átomos das diferentes camadas que o rodeiam. 


Fig.9 — Módulo |F (r)] da transformada complexa da função 
x (k). A tracejado está representada uma função “janela”, J(r). 


O primeiro máximo importante corresponde à contribui- 
ção da primeira camada em torno do átomo absorvente. 
Outros máximos de menor grandeza, acima do fundo, cor- 
respondem a camadas de vizinhos mais distantes. 

A análise de Fourier da função x (k) pode ser refinada 
da seguinte maneira: faz-se a transformada de Fourier in- 
versa para cada máximo da função de estrutura separada- 
mente. Para isto “filtra-se” a função de estrutura por meio 
duma função “janela” J (r) que dá realce aos valores numé- 
ricos da função no intervalo em que se situa esse máximo, 
multiplicando-os por um factor que só tem valores apreciá- 
veis nesse mesmo intervalo. A função assim obtida, no es- 
paço de k, tem que ser corrigida (5) por ter sido obtida por 
filtragem. O espectro finalmente obtido, x" (k) corres- 
ponde à contribuição da referida camada de átomos vizi- 
nhos para o espectro de EXAFS. De x (k) pode derivar-se 
a função amplitude. 


2r, 
E 271 ———— 
At)= bm] e 2H * (9 
k Tó 
J 
O. (rm k) = 2k r(k) + 28, (k) + O (k) (16) 


pelo processo descrito em (41). 


O valor da distância interatómica é dado então por 


"= o; (r;, k) — 26", (k) — 2 (k) (17) 
2k 

Esta expressão não dá para r; (k) um valor independente 
de k devido às incertezas experimentais, às incorrecções 
cometidas na atribuição teórica dos valores às diferenças de 
fase, ou na escolha dos valores de E. Em 8 Â podem 
obter-se valores de r com um erro de + 0,05 À. Um me- 
lhor ajuste de E, permite atingir um erro apenas de 
+ 0,01 À. Estes valores foram obtidos por Martens et al (33) 

para compostos de Cu. 
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A. As Diferenças de Fase 


Como vimos, para determinar as distâncias interatómi- 
cas interessa conhecer as diferenças de fase. 

As diferenças de fase podem obter-se por determinações 
empíricas ou por cálculos bascados em considerações tcóri- 
cas. No primeiro caso obtém-se apenas a diferença de fase 
total y (k). No segundo são calculadas separadamente a 
diferença de fase do átomo central, ô, (k), e a diferença de 
fase da onda rectrodifundida, O (k). 

Os valores obtidos teoricamente encontram-se tabelados 
por B.K. Teo e P.A. Lee (45) para grande parte dos elemen- 
tos do quadro periódico e é possivel obter por interpolação 


QUADRO 4 


Valores das diferenças de fase do átomo central ô, (Kk) 
e da onda rectrodifundida (2 (k) para diversos elementos 
(valores calculados para k = 5, 1967 A) 


ô (k) O (k) 
(rad) (rad) 
O] =7, 1700 D,I85A 
Br 2,6538 59716 
In 6,3803 6,6205 


valores das diferenças de fase para átomos que não figuram 
nas tabelas. No quadro 4 figuram os valores obtidos direc- 
tamente, ou por interpolação, para o oxigénio, bromo e 
indio. 

Na fig. 10 pode observar-se a diferença de fase introdu- 
zida pelo átomo central para o iodo, bromo, cloro e fluor, é 
na fig. 11 a diferença de fase correspondente à onda rectro- 
difundida para os mesmos elementos. 


10.0. 


dl kl [rad.) 


Fig.10 — Dependência relativamente a k da função de fase do 


átomo central (absorvente), ô; (k), para diversos elementos (I, Br, 
Cle F). 


10.0 4 
8.0 4 


6.0 


2.0 


0.0 


Da lkl trad) 


e 8.0 


Fig.11l — Dependência relativamente a k da função de fase rectro- 
difundida, & (k), para diversos elementos, (I, Br. Cle F). 


As determinações empíricas baseiam-se numa hipótese 
chamada de transferibilidade quimica. Esta hipótese é vá- 
lida porque para energias cinéticas suficientemente eleva- 
das, acima de 50 eV (onde o processo de EXAFS é domi- 
nado pelos electrões mais internos) as diferenças de fase são 
insensíveis às modificações da estrutura interatómica, As- 
sim, para um cristal, de referência, podemos extrair o valor 
de y (k) e utilizá-lo para o mesmo par de átomos num 
sistema em que Wy (k) é desconhecido (por exemplo para 
uma solução aquosa). A validade deste processo pode de- 
monstrar-se comparando o valor duma distância interató- 
mica conhecida com a obtida com base nesta hipótese. 

Separadas as contribuições da fase 2 .(k), relativa a uma 
substância conhecida, pelo processo descrito em (41) obt- 
êm-se 


Do (k) = 2k re + be (k) (18) 


Separando da mesma maneira a fase D 4 (k) para a subs- 
tância a investigar, e subtraindo membro a membro, ob- 
tém-se 


8d (k) — Be (k) == 2k (rg — Te) (19) 


Do coeficiente angular da recta obtém-se rg. 


B. A Amplitude da Onda Rectrodifundida 


A amplitude da onda rectrodifundida pode determinar- 
-se teórica ou experimentalmente, neste último caso fa- 
zendo intervir a hipótese de transferibilidade (entre um sis- 
tema conhecido e um desconhecido). 

Os valores calculados satisfazem qualitativamente, mas 
quantitativamente por vezes sobrestimam a amplitude. 

Dois tipos de circunstâncias da difusão inelástica ocor- 


rem quando se dá o fenómeno de EXAFS que conduzem a 
uma redução da amplitude relativamente ao valor calcu- 
lado teoricamente. As do primeiro tipo são excitações múl- 
tiplas no átomo central pois o excesso de energia E-E, no 
processo de fotoionização pode excitar outros electrões 
mais fracamente ligados. Em geral estas excitações múlti- 
plas têm um espectro largo que tende a desvanecer o espec- 
tro de EXAFS. As do segundo tipo são situações de falta de 
coerência da onda rectrodifundida pelos átomos vizinhos 
relativamente à onda que parte do átomo central por ser 
obviamente finita a vida média do estado excitado: uma vez 
que o fenómeno de EXAFS resulta da interferência junto 
do núcleo da onda de matéria associada ao fotoelectrão, 
que parte do átomo central, com a onda de matéria rectro- 
difundida pelos átomos vizinhos, se o estado excitado do 
átomo central declinar antes do tempo de trânsito, a função 
de onda perde coerência e não se dá interferência. O termo 


Fig.12 — Dependência relativamente a k da função amplitude rec- 
trodifundida, |f, (k, 7 )|, para diversos elementos: a. F, Cl, Bre I: 
b. Br, Ge, Cu, Fe, Cre Ti, 
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do livre percurso médio À na equação de x (k) dá conta 
deste efeito de tempo de vida, 

Ora as amplitudes da onda rectrodifundida tabeladas por 
B.K. Teo e P.A. Lee (45) têm em conta esta falta de coerên- 
cia entre as funções de onda mas não têm em conta a multi- 
plicidade de excitações, o efeito a muitos corpos, dos elec- 
trões mais exteriores do átomo central. 

Um processo de corrigir o valor da amplitude consiste 
em introduzir na expressão de x (k) umfactor S(k) que é 
aproximadamente igual a | para os baixos valores de k e 
varia entre 0,6 € 0,8 para valores de k < 7 Á'!, 


QUADRO 5 


Valores da amplitude da onda rectrodifundida | f, k, m) | 
para diversos elementos. (Valores calculados para 


k = 3.7795 À) 
Z f, (k, 7) 
Ó BR 0,704] 
Br as (1,3836 
In 49 1,0478 


Na fig. 12 estão representadas funções amplitude 
|f, (k, m)| para o bromo e para outros elementos. 


No quadro 5 podem ler-se os valores da amplitude da onda 
rectrodifundida para o oxigénio, o bromo e o índio tabelado 
por B.K. Teo e P.A. Lee que como já referimos incluem o 

ne 

factor exp ( — — 

A (k) 
sim desnecessário o cálculo de À pelo processo descrito 
por exemplo por Stern et al. (23) ou mais recentemente por 
B. Lengeler e P. Eisenberger (37), que propõem um método 
para determinar os factores de perdas inelásticas para cada 
camada e para cada par de átomos absorvente-rectrodifu- 
sor a partir de compostos de referência de estrutura conhe- 
cida. 


) da expressão de x (k). Torna-se as- 


C. O Número de Átomos Vizinhos 


A função amplitude obtida da transformada de Fourier 
inversa, “filtrada” para cada camada, contém informação 
sobre o número de vizinhos dessa camada em torno do 
átomo central. Para se obter o valor de N, admite-se uma 
hipótese de transferibilidade quimica para a amplitude, 
pelo que é necessário dispor dum sistema de referência em 
que M, seja conhecido, 

Da expressão (1), considerando a amplitude para uma só 
camada 


2r 
2a o 
—2Zke a À (k) 


EE TI] É o 


A (k) = (20) 


kr? 


e calculando o logaritmo da razão das amplitudes para os 
dois sistemas, obtém-se 


Ac (k) : à 
Ind —) =.In | O) + 2 (04 — 07) 
Au (k) Na fá y 
r r 
+ 2 ARES: =—=.3 (21) 


Admite-se que Ag (k) = Ac(k)jeque(ra -r)<<A (k), 
nestas condições o terceiro termo é despresável. Da repre- 
sentação gráfica do primeiro membro da expressão (21) em 
função dek?, pode obter-se o valor de Np e og depois de 
obtido ry por um dos processos descritos atrás. 

Ao contrário do que se passa com a transferibilidade das 
diferenças de fase, a transferibilidade da amplitude pode ser 
afectada por vários factos. Esses factos podem ser a má 
determinação empírica de A (k), se a gama de valores 
de k de que dispomos for pouco extensa, ou a circunstân- 
cia da relação sinal ruído poder ser pequena, ou ainda a 
incerteza na determinação dos factores de perdas inelásti- 
cas como foi referido em 5.B. 

O facto de a transformada de Fourier ser filtrada intro- 
duz uma distorção na amplitude. Uma correcção a esta 
distorção faz-se usando o valor teórico da amplitude para 
calcular x (k) que será analisado da mesma forma exacta- 
mente que a função x (k) obtida experimentalmente, e 
comparando as duas, variando da mesma forma os valores 
limites dos kk nas duas funções (a experimental e a teó- 
rica). 

Para a primeira camada de átomos vizinhos em torno do 
atomo central é possível no entanto obter valores de N, com 
um erro inferior a + 10%. 

A análise dum espectro de EXAFS dá-nos além dos va- 
lores das distâncias interatómicas, informação sobre o nú- 
mero de coordenação e ainda sobre um terceiro parâmetro, 
dificil de obter doutra forma, o coeficiente de Debye-Wal- 
ler, o, que é uma importante caracteristica da estrutura 
dum sistema, 


6. VANTAGENS DE EXAFS EM RELAÇÃO 
A OUTRAS TÉCNICAS DE INVESTIGAÇÃO 
DE ESTRUTURAS MOLECULARES 


As técnicas de difracção, dos raios X ou dos neutrões, 
aplicadas a substâncias que apresentam uma ordem ex- 
tensa, como os cristais, dão uma informação tridimensional 
da sua estrutura; porém, aplicadas a substâncias que têm 
apenas uma ordem local, como é o caso dos líquidos e dos 
sólidos amorfos, a informação estrutural é apenas a uma 
dimensão. Com efeito, a transformada de Fourier da inten- 
sidade de uma figura de difracção dá uma função de distri- 
buição radial a uma dimensão, que para cada distância dá 
todas as correlações que no espécimen ocorrem a essa dis- 
tância, sejam quais forem os pares de partículas que lhes 
dão lugar, enquanto que a técnica de EXAFS fornece uma 
informação selectiva, isto é, para cada espécie de átomo 
separadamente — desde que se observem sucessivamente as 
regiões espectrais em que se situam as descontinuidades de 
absorção dos diversos tipos de átomos do espécimen. Tal 
informação, o número e a espécie de átomos vizinhos e as 
suas distâncias ao átomo absorvente, está contida numa 
função de distribuição radial a uma dimensão, centrada no 
átomo absorvente. 

Embora a difracção dos raios X ou a ressonância magné- 
tica nuclear dêem maior precisão para determinar números 
de coordenação, a EXAFS permite obter na determinação 
das distâncias entre um átomo ou um ião e as suas primeiras 
camadas de coordenação um rigor de 0,01 À. O acordo entre 
os resultados (distâncias interatómicas) obtidas pela difrac- 
ção dos raios X e por meio de EXAFS é bom para os cristais 
e difere de 2% nos liquidos. À discrepância provém de ser 
diferente a região do espaço dos momentos k: na difracção 
dos raios X o parâmetro de difracção k varia entre 1 e 
18 À"! enquanto queem EXAFS6 À! < MW <2ÃÁ. 


Além do que atrás se disse, um aspecto notável em 
EXAFS é ser possível conhecer a espécie de átomo coorde- 
nado pelo átomo central (pelo átomo em torno do qual 
estamos a estudar a estrutura) e se a camada de coordena- 
ção for constituída por duas espécies de átomos diferentes, 
poder distingui-los. Isto é possível através do perfil da am- 
plitude da onda rectrodifundida (obtida por transformada 
de Fourier inversa, filtrada para uma camada) sobretudo se 
o número atómico Z dessas duas espécies de átomo for 
bastante diferente (A Z > 4). 

A espectroscopia de EXAFS é de facto uma técnica privi- 
legiada para o estudo da estrutura molecular de liquidos, de 
polímeros, de semicondutores, de substâncias não cristali- 
nas que interessam não só aos domínios da Fisica e da 
Química mas também à Biologia. 

Agradece-se à American Physical Society ter autorizado 
a reprodução das figuras 2, 3,4 e 5 da Physsical Review B, 
17 (1978) 1481 e das figuras 7, 8 e 9 da Reviews of Modern 
Physics, 53 (1981) 769. 

As figuras 10, 11, 12a e 12b foram reproduzidas com 
autorização da revista Journal of the American Chemical 
Society 101 (1979) 2815, copyright (1979) American Che- 
mical Society. 

Este artigo foi redigido com base no projecto de investi- 
gação que a autora apresentou em Setembro de 1983 à 
Faculdade de Ciências da Universidade de Lisboa para 
prova complementar do doutoramento em Física. O tema 
do projecto foi escolhido pela autora de entre alguns que o 
Professor M. Alves Marques lhe propôs. 


BIBLIOGRAFIA 


(1) International Critical Tables of Numerical Data, 
Physics, Chemistry and Technology, vol VI, pg. 33, 
ed. E.W. Washburn (Macgraw-Hill, 1928). 

(2) W.M. Dumond and P. Kirkpatrick, Rev. Sci. Instr. | 
(1930) 88 citado em “X-Rays in Theory and Experi- 
ment” A. Compton & S. K. Allison, pág. 750. 

(3) W.L. Bragg and J. West, Phil, Mag. IO (1930) 823 
citado por (5) 

(4) T. Johansson, Naturw. 20 (1932) 758. 

(5) J. Waser and V, Schomaker, Rev. Mod. Phys. 25 
(1953) 671. 

(6) F. Debot, Physica 21 (1955) 605. 

(7) F.W. Lytle, Phys. Non-Cryst. Sol. (1965) 12 ed. J.A. 
Prins (North Holland, Amsterdam). 

(8) F.W. Lytle, Advances in X-ray Analysis, 9 (1966) 398, 
ed. G.R. Mallet, M. Fay, and W.M. Mueller (Plenum, 
N-Y). 

(9) J.A. Bearden, Rev. Mod. Phys. 39 (1967) 78. 

(10) JA. Bearden and A.F. Burr, Rev. Mod. Phys. 39 
(1967) 125. 

(11) International Tables for X-ray Crystallography, ed. K. 
Lonsdale (Kynoch Press) 1968 vol III. 

(12) €. de Boor, J. Approx. Th. | (1968) 219 e 6 (1972) 50 
citado por (41). 

(13) J.B. Pendry, J. Phys. €. 2 (1969) 1215. 

(14) L. Salgueiro e J.G. Ferreira, “Introdução à Física 
Atómica e Nuclear”, vol I, pág. 197, (1970). 


(15) D.E. Savers, F.W. Lytle and EA. Stern, Advances in 
X-ray Analysis, 13 (1970) 248, ed. B.L. Henke, J.B. 
Newkirk and G.R. Mallet (Plenum, N-Y). 

(16) D.E. Savers, E.A. Stern and F.W. Lytle, Phys. Rev. 
Letters. 27 (1971) 1204. 

(17) D.E. Savers, F.W. Lytle and E.A. Stern, J. Non-Cryst. 
Sol. 8-10 (1972) 401. 

(18) E.A. Stern and D.E. Sayers, Phys. Rev. Letters 30 
(1973) 174. 

(19) E.A. Stern, Phys. Rev. B 10 (1974) 3027. 

(20) FW. Lytle, D.E. Sayers and E.A. Stern, Phys, Rev. 
B 11 (1975) 4825. 

(21) P.A. Lee and J.B. Pendry, Phys. Rev. B 11 (1975) 
2795. 

(22) C.A, Ashley and S. Doniach, Phys. Rev. B 11 (1975) 
1279. 

(23) E.A. Stern, D.E. Sayers and F.W, Lytle, Phys. Rev. 
B 11 (1975) 4836. 

(24) B.M. Kinkaid and P. Eisenberger, Phys. Rev. Letters 
34 (1975) 1361. 

(25) P. Eisenberger and B.M. Kinkaid, Chem. Phys. Let- 
ters 36 (1975) 134. 

(26) P.A. Fox, A.D. Hall and N.L. Schryer, Bell Laborato- 
ries Computing Science Technical Report n.( 47 (1976) 
citado por (41), 

(27) 6.8. Brown, P. Eisenberger and P. Schmidt, Solid 
State Com. 24 (1977) 201. 

(28) P.A. Lee and G, Beni, Phys. Rev. B 15 (1977) 2862. 

(29) D,R. Sandstrom and H.W. Dodgen, J.Chem. Phys. 67 
(1977) 473. 

(30) E.A. Stern, Cont. Phys. 19 (1978) 289. 

(31) P. Eisenberger and B.M. Kinkaid, Science 200 (1978) 
1441. 

(32) P.H. Citrin, P. Eisenberger and R.€. Hewitt, Phys. 
Rev. Letters 41 (1978) 309. 

(33) G. Martens, P. Rabe, N. Schwentner and A. Werner, 
Phys. Rev. B 17 (1978) 1481. 

(34) A. Fontaine, P. Lagarde and D. Raoux, Phys. Rev. 


Letters, 41 (1978) 504. 
(35) B.K. Teo and P.A. Lee, Am. Chem. Soc. 101 (1979) 
2815. 


(36) D,R. Sandstrom, J. Chem. Phys. 71 (1979) 2381. 

(37) B. Lengeler and P. Eisenberger, Phys. Rev. B 21 
(1980) 4507. 

(38) E.A. Stern, B.A. Bunker and S.M. Heald, Phys. Rev. 


B 21 (1980) 5521. " 
(49) P. Eisenberger and B. Lengeler, Phys. Rev. B 22 
(1980) 3551. 


(40) P. Lagarde, A. Fontaine, D, Raoux, A. Sadoc and P, 
Migliardo, J, Chem. Phys. 72 (1980) 3061. 

(41) P.A. Lee, P.H. Citrin, P. Eisenberger and B.M. 
Kinkaid, Rev. Mod. Phys. 53 (1981) 769. 

(42) W. Thulke, R. Haensel, and P. Rabe, em “EXAFS 
and Near Edge Structure” — Proceedings of the In- 
ternational Conference, Frascati — Italy, 13-17 de Se- 
tembro 1982 (Springer — Verlag, 1983). 

(43) G. Licheri, G. Pinna and G. Navarra Z. Naturf, 38 
(1983) 559. 

(44) G. Licheri, G. Paschina, G, Piccaluga and G. Pinna, J. 
Chem. Phys. 79 (1983) 2168. 


87 


A GENERALIZED INTEGRAL CHARGE RELATION 


FOR BIPOLAR DEVICES 


P.L. BORGES TEIXEIRA 


Centro de Electrónica Aplicada (INIC) e Dept. Eng. Elect. Instituto Superior Técnico 


RESUMO 


Apresenta-se uma generalização original da relação 
de Gummel, que faz intervir a condução ambipolar em 
todas as regiões do dispositivo. Com base nesta relação, 
propõe-se um modelo de carga integral generalizado, 
que se compara com as relações de Gummel e Poon. 
Obtêm-se conclusões teóricas sobre a injecção em alto 
nível, e sobre a simetria da corrente de transporte. 


ABSTRACT 


A generalized, new version of the Gummel relation 
is presented. Ambipolar conduction is assumed in all 
transistor regions. A generalized integral charge model 
is proposed, and compared with the Gummel and Poon 
relations. Theoretical conclusions on the high-level 
injection and symmetry of transport current are deri- 
ved. 


INTRODUCTION 


In a classical paper [1], Gummel presented an 
integral charge relation which shortly became a valua- 
ble tool for the one-dimensional modelling of bipolar 
devices [2], [3]. However, the low-level approxima- 
tions made in deriving that relation do not generally 
apply to normal, high-level operation of bipolar power 
transistors. The aim of this paper is to briefly describe 
a new generalized integral charge relation where no 
such assumptions are introduced ab initio. 


BASE REGION 


The problem is restricted to one-dimensional, x 
directed planar current flow under static conditions. 
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We shall assume the following approximations which 
by no means are restrictive to low-level operation: 
macroscopic formulation of transport equations; no 
degenerate concentrations, 1. e. the Einstein-Nernst rela- 
tion is valid, and bandgap narrowing effects are neglec- 
ted; contact surfaces are in equilibrium; quasi-neutra- 
lity of intermediate regions. Considering an ambipolar 
current flow, we may derive from the basic transport 
equations the following general relation 


| PO) nO). Xo> " "xo Jnlx) so d 
na X, =) X, q D,(x) n; dias ad 
-[. x) d | 
“ ED n(x) dx (1) 


where the symbols in the above equation have their 
usual meaning, and x; and x, are any two values of x 
inside the semiconductor structure. We shall first apply 
this relation to the base quasi-neutral region (0,Ws) of 
a NPN transistor. Introducing the charge densities 


We - 
dons fo “q pGo) dx 

Ro | 
a f 0 “q n(x) dx (2) 


and making use of the continuity equations, we obtain, 
after some involved calculations, for the electron cur- 
rent at the BC boundary plane 


(pn/n))ws ii (pn/n9, 
ERR ra 


J(WeB) = Bans) dz 
+ (Wo) b ER + Ro(Ws) (3) 
Ta 


Here Re(Wes) stands for an ambipolar recombination 
current, and b = &,/&p. Is is also possible to derive 
the electron current at the EB boundary, which is 
given by 


(pn/nwa — (pn/nd, 
ia 


AD) = pe eg 
Qps 


+ Jo(O) b E + Re(0) (4) 
p 


On 
Equations (3) and (4) are integral charge relations 
that describe the ambipolar current flow through the 
base. Hole currents Jo(Ws) and Jp(0) are injection com- 
ponents respectively into the collector and the emitter. 
They both express boundary conditions which depend 
on the charge transport on those regions. It can be 
shown that recombination currents Rg are composed of 
two terms. One is essentially independent of the injec- 
tion level, and the other, proportional to (bQ,s/Qps) 
is only important under high-level conditions. The 
component Ja(Ws) - (bQne/Qpe) in eg. (3), or its coun- 
terpart in eg. (4), has the physical meaning of an 
electron (minority) drift current due to the electric 
ficld produced by the hole (majority) current. The 
transport currents retain their usual signification. 


EMITTER REGION 


This region extends from x = — We (contact) to 
x = — Wº (EB junction). Integrating (1) between these 
limits, we obtain after some manipulations 


Em. (pn/ni) w —1] 
es Jo( = We )=q'Dpenir o | 


, À dee 
Es Medos + Re (5) 
where 
— WE — WE 
Que af qn(x) dx ; QF + [ q | nn | dx: 
— WE — WE 
— Wi 
QpE a | qp(x) dx (6) 
cons 


COLLECTOR REGION 


This region extends from x=Wí.to x=W. defi- 
ned by the equilibrium condition p(Wo) . n(We) = 


= n;. The dependency Wc(Jc,Vce) is mainly required 
at the quasi-saturation. We shall assume Wc as cons- 
tant. Integrating (1) between the above limits, we 
obtain 


(pn/nidwe — 1 


EWE) = er re q Dent + 
Que E 
+ Dão We j= Re (7) 
where 
We Wc 
Gu a [ anqudx ; & a ) q | n6)-n0o |ax; 
c c 
Wo 
Qpc 4 f gp(x)dx (8) 
W' 
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GENERALIZED CHARGE-CONTROL EQUATIONS 


In order to link the current boundary conditions 
(5), (4), (5) and (7), the transition regions across EB 
and BC junctions are taken as in quasi-equilibrium, 
and recombination neglected. We may thus write 


it WO=T, O 5 Ti WO=T, (W,) (9) 


Eliminating Jo(Ws) and Jp(0) from the above men- 
tioned system, the x directed emitter current Je. = — 
— |r(— We) and collector Je=Ir(Wc) current are 
given by 


l 
I—ô%gôs 


1 ph 
+ (1 + 3p)Je (5) dk 
ço 


+ (1+)m (2) 
(o, 


1 WB 


JEx — 


( — [ + d)Jei + 


+ Tres (10) 


| — 8côp 


+ (1 + ôs)Jci (E) e 
n. WEB 


Jcx = e] K + 8c)Ja + 


pn 
— (1 + ôc)Ja: (7) +) (11) 
n RCx 


RR 


Here Je, Ja, and Jo are reference currents, J'Rex 
and J'Recx indicate a condensation of recombination 
terms, and the coefficients 
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á bQne i QpB i bQnc 


provide an indication of the high-level injection into 
the regions they refer. Let us consider the (1 — 87ô5)' 
factor cf Jex, and suppose that ô; and ôs may simulta- 
neously assume finite, non-zero values. The result 
d:ôB = | may then well take place, and thus Jp”. 
This is absurd, and reveals an inconsistency at equation 
(10). We thereby conclude that high-level injection (as 
described by ôE, às, Or dc) can only occur into the region 
with lower impurity concentration. For the forward- 
-biased EB junction, we thus have dr <ôs. The above 
conclusion may be extended to the (1 — 3,8) ! factor 
of Jcx. 

To make use of (10) and (11) on device modelling, 
it is convenient to normalize charges Q pm to some refe- 
rence of equilibrium values (Qp,n)o. The normalized, 
dimensionless charges may then be written as 


QB O pE PC 


Que =— "= + Gp — 
QB, QuE, 


(13) 


Considering the neutrality relations for each quasi- 
neutral region, we obtain the following generalized 
mitegral charge equations 


Jes | b+1 
1 + 


Jes | a” b+1 / pn ) 
ii GE : b me nº WB 
l b+1 


a A (1 + b QpE ERR qu | Jres (14) 


Jas b+1 
Jet) = qe |( 1 + — que ) + 


Jos [ooo pn 
Cuemsnse (E) = 
Jas | n / wB 
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| b-+1 | 
de amis ( + qe) I+(b+1) dos | Jrex (15) 
é * 


Que — L + Qpe + Qua (16) 
qee = | + que + Que 


and Jas, Jes, and Jces are regional saturation currents. 
As in the Gummel and Poon (GP) model, we may 
define J,./Jas = Bp' and Joc/Tas = B;, etc. A compa- 
rison of equations (14) — (15) with the Gummel rela- 
ton shows that the integral charge is now replaced 
by composite integral charges qsr or qec. This corres- 
ponds to the logical extension of Gummel's integral 
charge concept. As a direct consequence of (14) — (15), 
the base current Js=Jcx—Jex can now be analitically 
defined, in contrast with the GP model, where its 
introduction is mainly empirical. Another interesting 
feature is that the transport current is no longer sym- 
metric, 1. e., not common to the emitter and to the 
collector currents, unless que=qpc=0, which means 
low-level injection. But in this case, equations (14) — 
— (15) reduce to the Gummel relation, where only 
charge transport through the base is considered, and 
taken as symmetric. In a planar power device, where 
Queo> QproP Quco, We may assume that q: =0; however, 
Gpc>1 inside the saturation or quasi-saturation ope- 
ration. 

The generalized integral charge equations may be 
used at the modelling of bipolar devices, such as tran- 
sistors, thyristors and pn or psn diodes. A model for 
computer simulation of high-power transistors was 
develloped, and some related results were presented 
elsewhere [4]. 
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COORDINATES TRANSFORMATION FROM THE 
PRINCIPLE OF RELATIVITY 
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SUMÁRIO 


Apresentamos neste trabalho uma dedução das leis 
de transformação das coordenadas cartesianas ortogo- 
nais associadas a referenciais de inércia. Apenas se 
recorre ao princípio de relatividade. Assim se obtém 
uma forma geral dependente de uma constante, “a” 
determinável a partir do conhecimento da velocidade 
vems (xi) e sua transformada v em S (x). 

A obtenção de “a” é, portanto, matéria de expe- 
rência directa ou de hipótese consistente com anterio- 
res factos experimentais. 


SUMMARY 


We present a derivation of the transformation law 
of ortogonal Cartesian coordinates associated with 
inertial frames. Using only the principle of relativity, 
we have obtained general transformations depending 
on a constant “a”. The determination of “a” is possi- 
ble after the knowledge of a velocity vin S (x') and the 
transformed value v in S (x!). So, the correct value of 
“a” is a matter of a direct experiment or results from 
a hypothesis consistent with experimental facts. 


The space-time transformations of Special Relativity 
result necessarily from the two basic assumptions: 


| — Equivalence of inertial frames in uniform rela- 
tive motion: special principle of relativity; 

2 — Independence of light velocity from source 
inertial frame and observer inertial frame: in 
the Special Relativity this velocity is a univer- 


Pr PF 


sal constant Cc. 
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From these two postulates we may conceive 
different situations from which result the Lorentz 
Transformations [1]. 

If we use only the first postulate it is possible to 
obtain a more general form of space-time transforma- 
tions including Lorentz transformations or any other. 
The correct option is a matter of experience. 

This fact has been pointed out by others, namely 
Aharoni [2] and Lee-Kalotas [3]. 

The aim of this paper is the same and we believe 
that the new fact is the derivation of the correct law 
on the basis of an experimental determination of a 
constant, depending on the knowledge of the two cor- 
respondent values of the velocities [v, S (x)] and [v, 
S (x)]. We have obtained the analytical expression of 
this constant. On the other hand, we think, the deriva- 
tion of the results is particulary simple. 


The homogeneity and isotropy of space-time are 
implicitly contained in the principle of relativity res- 
tricted to inertial frames [4]. 

From the homogeneity — translation invariance — 
we can impose the same coordinates to the origins of 
the space-time: O [0,0,0,0] = O [0,0,0,0]. 

From the isotropy — rotation | invariance — we 
can impose at t=o, the coincidence of the homonymous 
axis of ortogonal Cartesian coordinates: 


t=o0, XXxEXx, YyEyy, Z2 E Zz. 


We supose now a uniform relative motion paralle) 
to the x-axis: relative velocity O/S = v, ex. 

Under the circumstances we are led, [2], [4], to 
the set of linear transformation equations: 
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x=ax+ bt 

t=ecx+dt (1) 
y=y 

Z=2Z 


So, the crucial point is the determination of [a, b, 
c, dj. 
From (1) we obtain 


b ad—bc - 
== dm 2 
x d t + d X (2a) 
aus o l 
ANNE + (2b) 
a a 


By the principle of relativity the character of the 
motion observed is relative, so 


Velocity O/S = — Velocity O/S 
(3) 


The coordinates of O in S are S (0,0,0:t) and the 
coordinates of O in S are S (0,0,0,t). Then, if we take 
in (2a) x = o we have x/t = vr. 

Similarly, if we take x = o in (2b) we shall have 
X=—w. 

Consequently we obtain the following results 


a o a=d (4) 
d a 
and the equations (2) take the form 
2 —- bc — 
x=wtt————x (5a) 
- - 1 
A MR (5b) 


Suppose now a bar of length x in x-axis of S. By 
(5a) and at the time t=o, we have in referential frame 
5 a correspondent length x. These two lengths are in 
ratio 

x a — bc 
x a 


The same experiment for a bar of length x in S 
leads in S and at the time t = o to a correspondent 


value x. We have now 


+ (6b) 


[a 
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By the principle of relativity this ratio must be 
exactly the same since the experiments are conceptually 
identical and made for an equal modulus of the relative 
velocity, vr. 

Consequently we must have 


a — bc 2. É 


=—— +g-te=] (7) 
a 
and the equations (5) become 
| E = 
x= vt+ o X (8a) 
- - 
à cm (8b) 
x 
* x 


These are the more general transformations we 
may write using only the first postulate: the principle 
of relativity limited to inertial frames. 

The determination of “a” is a matter of experience. 


* 


In a general way we may conclude from (1) that 
x 
E 


for = y = const. there corresponds a unique 
X 


v = const. which means that a point with 
É ia O TES aa cui e. 

coordinates P(x, y, z, t) moving with a velocity v = — 
a L 

= const. in Sis seenin S as a point P (x,y,z,t) moving 


X 
with velocity v = — = const. Obviously we have 
— L 


y=yez=z 

The movement of P may be simulated or may cor- 
respond to a natural movement of a particle or group 
of waves. 

In any way it is a conceptual and possible experi- 
ment with different physical means. 
If we know a pair v=>v this is sufficient to find “a”. 


If we take in (8a) and (8b), x= vt and x = vt 


we obtain 


= wt=9" (9) 

” - 1 

(v+ v)t= i vt (9b) 
from which results 

+ Fe É (10) 


(v — vo) (v + vo) 


So and if the principle of relativity is observed in 
the physical world, accurate measurements of v must 
show the existence of a constant “a” depending only on 
the relative velocity vr. 

Once more by the principle of relativity “a” must 
he an even function of v.: a(v;) = a(— vi). 


* 


The equations (8a) and (8b) may be written in the 
equivalent forms: 


x=a(x— vt) x=a(x+ yt) 
e É a: | a E 
isa(t-*D q) t=af+t—— d 
a Vr a Nr 
(lla) (11b) 


By the principle of relativity we must obtain (11b) 
from (1la) with a simple change of the signal of v; 
which implies in fact a(v;) = a(— vo). 


* 


The Mechanics of Newton don't put any limitation 
on a velocity of a material particle. In principle and 
in the scope of this theory we may accelerate a particle 
uD to vy = oo. On the basis that v > v, for v; > 0, we 
obtain too v = co and from (10) there result a” = 1 
and à = *l, 

Since we must have the identity transformation 
for v. = o we adopt a = 1. 

Introducing this value in (11b) we obtain the Gali- 
lilean transformation 


(13) 


V=v+v 


On the basis of the Fizeau's experiment in a mo- 
ving stream of water, [5], we have: 


— c 


c 
-—-=——. a v=——- + (1 -— 
n n né 


where c is the velocity of the light in vacuum and n 
the refractive index of the water; dispersion effects are 
neglected. 


) Vr 


If we introduce v and vy in (10) we obtain 


Ae-p-keB-R E (14) 
a 


where 8 = v;/c and k = (n — =), 


Neglecting the higher order terms in 8 we obtain 


2 
Yr 


- (15) 


l 


This is an interesting example as by experience we 
may reveal that a 1 which is an essential fact be- 
cause a = 1 cancels the desyncronization and the time 
becomes absolute, t = t. 


This is the right moment to appeal to the second 
postulate of the Special Relativity. 
From this we obtain a convenient pair for our 


purpose: v=v=c. 


Introducing in (10) we obtain 


l | 1 
as É e DE E 16 
aí l / y ( ) 


that corresponds to the Lorentz transformation 


x=y (X+vw t) 


(17) 
a WE a 
t=y 6+5-») 
C 
and an associate law of velocity transformation 
v+v | 
v= AR A (18) 
| + vw/C 
de 
* * 


Finally we note that it is always possible to simulate 
in S a non-intrinsic causal process where the velocity 
ot Pis v = co, For example, this can be done with 
a line of lamps triggered by a wave front parallel to 
this line. 


But if take v = oo in (18) we obtain a finite value 
2 
c 
V=— >e. 
Vr 
So we are not able to realize an artificial process 
where on both frames we have v= v = 00. We be- 
lieve this is an essential point. 
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Conclusions: 


— We have established a general form of coordi- 
nates transformation using only the principle of 
relativity and concerning ortogonal Cartesian 
coordinates. 

In that form we have a degree of liberty: it is 


vv 
the constant &* = ———— 


(v— ve) (v+vr) 

— It is sufficient to know two velocities, v and 
v, to determine “a”. This is a matter of expe- 
rience. 

— The second postulate of the Special Relativity 
provides a convenient pair v = v=c leading 
to a = 1/(1 — 8?) which corresponds to the 
Lorentz Transformations. 
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— We can use another kind of phenomena to deter- 


Fr Pr 


mine “a. 
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EXPERIMENTAL RESULTS OF THE CONDITIONAL 
SIMULATION OF AN IRON OREBODY 


ANTÓNIO JORGE DE SOUSA 


ABSTRACT 


A numerical model ofa zone in an iron orebody is obtai- 
ned by the conditional simulation of Fe grades in technolo- 
gical blocks. This technique ensures that simulated values 
match the mean structural features of the mineralization 
(auto-correlation and histogram) and assume experimental 
grades in sampled points. 


RESUMO 


Através da simulação condicional dos teores nos blocos 
tecnológicos de uma zona de um jazigo de Ferro, obtém-se 
um modelo numérico que reproduz os aspectos estruturais 
da mineralização (auto-correlação e histograma) e toma os 
valores reais nos pontos amostrados. 


In the planning stage of a minig project, it is often useful 
to establish a model of the orebody, where the variability of 
a quality parameter of the ore is reproduced. The fluctua- 
tions of this parameter in a certain time basis may then be 
predicted from the model, allowing the assessment of the 
impact of these flutuations on the mill feed so as to define 
blending policies (if necessary). The mining method opera- 
tor may also be applied to an input provided by the model, 
im order to study the effect of extracting procedures on the 
quality parameter variability and compare different exploi- 
tation alternatives. 

This model of the orebody, which must provide unbiased 
variances, is constructed by the Conditional Simulation 
technique (JOURNEL, 1974). If the sampled values of the 
quality parameter are viewed as the realization of a Ran- 
dom Function (MATHERON, 1970), the Conditional Si- 
mulation technique generates another realization of the 
same Random Function, which reflects the structural fea- 
tures revealed by experimental data (auto-correlation and 
histograms and assumes the observed values in the sampled 
points (it is conditional to reality). 


Original recebido para publicação em 15/1/81. 


Theory of conditional simulation 15 presented elsewhere, 
MATHERON (1972), JOURNEL (1974), DOWD (1979) 
and SOUSA (1979). In this note, an aplication to Fe grades 
is performed, in order to establish a model of the variability 
of that parameter in technological blocks of an iron ore- 
body. 


EXPERIMENTAL RESULTS OF SIMULATION 


The Moncorvo iron orebody (Northest of Portugal) was 
recognized by a set of drillings, sampled every two meters 
for Fe grade. This is the parameter to be simulated and its 
spatial structure was depicted by the variogram, an 
auto--correlation function, definned by: 


y(h)= VAR [Z(x+h) — Z(x)] 


where Z (x +h) and Z (x) are values of grades in two points 
linked by a vector h. Coordinates of points (in Rº) are 
denoted x. 

The experimental variogram obtained in the deposit is 
shown in Fig. 1. An isotropic model with nugget effect and 
two spheric type nested schemes was fitted to the experi- 
mental curve. The simulation technique ensures that the 
new realization of the random function obeys to this vario- 
gram model, the parameters of whichare €C, = |8 (Go)? ( nug- 
get effect), C, = 15 (%)? (sill of the first scheme); C,-8 
(%)* (sill of the second scheme); aj=7m (range of the 
first scheme); a, = 29 m (range of the second scheme). 


The experimental histogram of Fe grades in samples is 
shown in Fig. 2. Also, point values of the variable obtained 
from simulation are expected to follow this histogram. 


The simulation was performed in a 200x200x40 m zone 
of the orebody (where production is planned do start). In 
this zone, real data to which simulation is conditioned, are 
80 samples located in 4 drillholes. The general purpose of 
this simulation is to model the variability of grades in a 
daily base block of 25x12.5x 10 m, in order to assess the 
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FIGURA 1 
Experimental 
variogram 

and fitted model. 


— VYARIOGRAMA EXPERIMENTAL 


— MODELO 


FIGURA 2 
Histogram 

of Fe grade 
in samples. 


10 


effect of the mining method and blending policies into the 
luctuations expected in the feed of the dressing plant. 


The conditional simulation algorithm applied in this 
case, includes the following steps: 


|. Inverse Gaussian transform of data values Z(x;) in 
order to obtain a normal distribution of a new varia- 
ble Y(x;). 

+. Computation of variogram y« (h) for variable Y (x,) 
and parameter estimation of the theoretical scheme 
adjusted; determination of the spheric type covariance 
tunction C(h)by C(h) = C(0) Yw th), where C(h) 
15 the 3 dimensional covariance of Y(x) and C(0) is 
the à priori variance of data, Simulated normal points 
in Rº are expected to follow this covariance scheme. 

3. Generation of 119 random points im a line with a 
| dimensional covariance €, (u), the parameters of 
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which are deduced from the 3- dimensional model 
C(h): €C, (u) is relared to C(h) by C, (u) = 
=0h dh Ch). 


Integration, in a regular mesh of 23 040 points, of 
particular contributions of cach one of the | —dimen- 
sional above mentioned simulations, performed in the 
lines joinning the mid-points of opposite edges of a 
regular icosahedron (turning bands method, in G. 
Matheron's terminology (MATHERON, 1972) ). 
Steps |. and 2. are necessary because this procedure 
produces, by the central limit theorem, normal distri- 
buted points. The output of this step 1s a set of non- 
-conditional simulated point values Yç (x). 

Direct Gaussian transform of simulated point values 
in Rº. in order to match the histogram of samples, 


- Clustering of 45 simulated points in each 25x 12.5x 10 m 


block and averaging grade values per block. The result 


of this procedure is a set non-conditional simulated 
average grades per block (Zc(V,)). 

7. Inverse Gaussian transform of non-conditional simu- 
lated values Zç (Vi) im order to obtain a normal dis- 
tribution of Yo (Vo). 

8. Calculation of a linear estimator (by Kriging, MATHE- 
RON, 1970), of deviations (from Y(x)) — Yç (x))), 
in all simulated blocks and summing these values to 
Ys(V,). The output of this step is a set of conditional 
normal values Ye (V.). 

9. Direct Gaussian transform of Ye (V,) in order to ob- 
tain Ze (V,). 


By means of the conditional simulation technique, a nu- 
merical model of grades in blocks of the above mentioned 
zone was obtained. This model is the set of values Zc(V,). 

Taking this model as an input, a variogram of simulated 
grades in blocks was computed. This variogram is shown in 
Fig. 3 and shows a good agreement to the experimental 
one, obtained from samples by a regularization procedure 
(MATHERON, 1970). 

In order to depict the distribution of simulated grades 
per block, an histogram of these values was computed (Fig. 
4). This histogram is the basis for further planning. 

The exploitation method and blending procedures may 
be simulated from the numerical model containing Ze(V,) 
values. 
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FIGURA 3 — Variagram of simulated grades in blocks. 
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FIGURA 4 — Histogram of simulated grades in blocks, 
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MODELLING OF AN INDUCTION MOTOR SPEED 
CONTROL WITH SLIP-POWER RECOVERY 
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Centro de Automática da Universidade Técnica de Lisboa (INIC) 
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ABSTRACT 


The sliding — contacts equation is used to develop 
a dynamic model for the induction motor with static 
inverter on the rotor circuit. 

Good agreement is obtained with experimental 
results m the steady-state. 

With that model, equivalent Park machine equa- 
tions are developed. 


O. INTRODUCTION 


A speed regulating scheme using a wound rotor 
induction motor and a static synchronous inverter pro- 
ves to achieve good results as far as efficiency and 
rehability are concerned, 

Fig. 1 represents the scheme analysed in the paper. 
The asynchonous machine, fed from the a.c. source, 


has its secondary circuit connected to an a,c./a.c. con- 
verter which is composed by a rectifier bridge and a 
synchronous inverter. 


In spite of many papers which appeared during the 
sixties and in which this system was analysed [1], 
[2], [3], a satisfactory modelling has yet to be found. 


In a recent paper dealing with this subject [4] a new 
approach was made to solve the above mentioned pro- 
blem, nevertheless it is the authors" opinion that the 
model presents some limitations. In this paper the sli- 
ding-contacts equation [5] applied to the rectified 
generator [6] allows the development of a dynamic 
model. The steady-state behaviour of the system is 
studied and compared to experimental results. 

Making a change of variables, the Park equations 
can be obtained. The differences between the repre- 
sentation of this paper and the one proposed in refe- 
rence [4] thereby become evident. 


la 
e Re L 


Fig. 1 — Inverter speed-control system 


Original recebido para publicação em 20/2/84. 
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|. MODELLING 


The stator behaviour can be represented by two 
quadrature windings «, 8, when the polyphase power 
supply is symmetrical [7]. The rotor and the rectifier 
bridge can be seen as an armature of a D.C, machine. 
The fictitious brushes of the commutation circuit are 
not in a fixed position, but they make an angle p that 
varies according to the load. 


Fig. 2 — Equivalent machine 


The «, B circuits are filiform and their equations 
derive from Faraday law. For the “winding a” it is 
necessary to apply an extension of the Faraday law to 
the commutation circuits [5]. 


So the system equations are 


” dy, 
Ui= 
iz [rÀ I, LE dt (1) 
RR | 
Up = Rig + dt (2) 
ER Ore : 
Us = Ra la + dt E . Ex (3) 
p =p (state variables of the system) (4) 
do, a 
E = n= ? De, 5 
L=] di SE (5) 
Neglecting the saturation effects we have: 
W cem W = ES L 2 RE 12 
E mo. 2 a + 2 Lp ty + 
l 2 . ! 
E 2 L, L, + M,, L, I, + M,, à, l, 
where 
Mo =Ma=—McosP Mo=Ma=-—Msin? 


The equations (1) to (3) are then: 


di, di, 
u, =R. ). +L—=—McosP 

“ dt 

dis di, . 
=R,1,+L— Edi = Mein — PM, Gone (7) 

di, di, dis 
UWu=RiktHli——Mcosf-—=—M sin Pp — 

dt dt Ro o 

+9,Mi sine — O, Mi, cos (8) 


The linear and sinusoidal hypothesis are not neces- 
sary for the sliding contacts equation but they simplify 
the model without significant errors. 

It is necessary to complement these equations with 
another one that enables us to determine the brushes 
position, which is defined by diodes in conduction. 
So, we associate the equivalent brushes position to the 
phases where the voltage is higher [6]. 

So, using a prove phase with voltage U' and posi- 
tion £”, we have: 


- 9UÍ 
( | ú E i 
dE =p 


UT 


| di, = dis 
= M e mm AR 
O M sin P di cos dt + 


+(P— 69) Lis +O M i, cos E + 


+6,M is sin € (9) 


The equations (5) to (9) form a mathematical model 
for the system, but it presents a disadvantage which 
results from the fact that the stator variables are time- 
-dependent in steady-state. To obviate this problem, we 
make a change of variables in the stator: 


[À = cos O. sin O, d 


B sn O,—cos O, q 


je 


Fig. 3 — Change of variables 
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With this transformation, the system equations 
becomes: 


dias = : 
Uds — R. lds + Lar e a, Ls las — 


di, 


—M cos 5 + (0, +48) M is sin 3 (10) 


dias ” ; 
Uas — R. RA + La Es 9, L; las = 


di, E ê | 
+ M sin Ta + (0.48) M à cos 8 (11) 


, dia dias 
um=R ihi+-+la —M cos é + 


dt dt 
dias 


a Osr M las sin 0 — 
É 


+ M sin 8 


— O M ig cos 8 (12) 


0 M si 5 dias M 5 dias 
= sin ER + cos dt + 


(9, + ô) L, Ee Bis M ia cos 8 + 
+ O M las sin 8 (13) 
M ja (ias sin 8 + ias cos 8) — 


do, 


raio O, (14) 


— [=] 


where 


0 =p— O, and O = 0,— 6, 

We choose the dg synchronous reference such that 
Ugs = O. This simplification is already classical when 
the induction machine is fed from a voltage source [8]. 

The voltage U, differs ffom U, (as we can see in 
fig. 1) by the presence of the R, L filter; and so, 


dis 


Up = Ry - 
id dt 


la + Le + U, 


The equation (12) is modified: 


P= Kr O e cos df 


da 


pi Ou M lgs cos ô (12.8) 
where 
Rr = Ra + Rr 
Lr = La + Lr 
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2. STEADY-STATE 


In this situation the currents in the equations (10) 
to (14) are constant and the equivalent brushes are 
fixed in the synchronous frame: 


So, the steady-state equations are given by: 
Us = R$ lu O, Ei 1a 4 6.M Lsin 8 (15) 
U= RB la=O lo M Loss (16) 


Up=Rr I.-0.M ly sin ô— 


0) — Os La I, Ene Da M Tas cos õ + 
+69. M IL sin ô (18) 


M I, (Iss sin ô + Ia; cos) — T,=0 (19) 


The resolution of the equations (15) to (18) is noí 
easy even in steady-state, so we solve this problem by 
a numerical method (Newton-Raphson). 

On the pictures that follow we represent the theo- 
retical and experimental results for a system defi- 
ned by: 


Us = V2.300V, R$= 32, Rr= 108, L.= 
=.46H, LL=.488H, M= 


=.44H,=50 Hz 
(1) Up=—-92V 
(2) Up = — 101.2 
(3) Up = — 193 
(4) Up = — 285 


Fig. 4 — Steady-state system characteristics 


In the above results we achieve a fair coincidence, 
so it not necessary make any correction [9], at least as 
far as the verified zones are concerned. 


3%. PARK EQUATIONS 


In the previously developed study, the rotor repre- 
sentation is obtained by means of a coil that takes a 
certain position in the dg chosen referencial. 

If we change the variables as follows, the m.m.f. 
of the rotor is generated by two circuits positioned on 
the same axis where the stator values are defined: 


jar = à cos 0 


q 


Fig. 5 — Change of variables 


According to this change, the equations (10) to (13) 
will be as follows: 


dias diar 
U 4 — R. ] 4 Ké EO de= CR 
) NESTE O a 
ll Mi (21) 
dias diar 
O=R ist =M—" — 
ias SF da ra aii 
Ed Õ, Eu lds + O. M ldr (22) 
U r— Ra Ê r Lo Ran Es Sigo dO TES 
; god in? Ra 
js Ou M laí + Os La lar (25) 
diar dias 
U C— | ra La =Ê GA 
e os laih dt dt 
pe Os M las 7: Or La ldr (24) 


Where 


Ui-= U. 2088 and Ur = — U. sin 5 


In Park equations the voltage U, is made up of two 
components Ua- and Usar whose values depend on the 
currents ia and igr. 


Since we have 


tg ô = — lar/ldr 


the above referred components are given by: 


Uar = U. cos (arctg [— iar/iar]) 
Ur = — U, sin (arctg [— iar/iar]) 


So we come to the conclusion that the component 
Ua: is not always necessarily equal to zero, as was 
done in the above mentioned paper [4]. 


4. CONCLUSION 


Our study enables us to conclude that the use of 
the sliding-contacts equation leads to an acceptable 
model to represent this system, at least in steady-state. 

The sliding-contacts equation apart from being 
easier to apply, has also proved to be power-ful means 
of analysis of the electromechanical conversion system, 
specially when compared to the classical method that 
results from the use of the Park transformation. 

We will try to develop this model in the future by 
presenting the dynamical behaviour of the system. 


REFERENCES 


[| — A. LAVI and R. 1. POLGE — «Indution Motor Speed 
Control with Static Inverter in the Rotor»; I.E. E.E. 
Transactions on Power Apparatus and Systems, vol. 
PAS-85, nº 1, January 1966. 

2 — WILLIAM SHEPHERD and JACK STANWAY — «Slip 
Power Recovery in an Induction Motor by the Use of a 
Thyristor Inverter»; I. E. E. E. Transactions of Industry 
and General Applications, vol. IGA-5, n.º 1, January/ 
!February 1969, 

3 — DEREK A. PAICE — «Speed Control of Large Induction 
Motors by Thrystor Converters»; 1. E. E. E. Transactions 
on Industry and General Applications, vol. IGA-5, n.º 5, 
September /October 1969. 

4— V.N. MITTLE, K. VENKATESAN and 5. C. GUPTA — 
«Stability Analysis of a Constant Torque Static Slip- 
-Power-Recovery Drive»; I. E. E. E. Transactions on 
Industry Applications, vol. IA-16, n.º 1, January/February 
1980. 

5—M. S. GARRIDO — «Les équations générales des machi- 
nes déduites de Velectromagnetisme»; Revue E, vol VI, 
n.º «9, 1971. 

6-—H,. BUYSE and M. S. GARRIDO — «A new dynamical 
model of rectificd output A. C. generator»; Electric Mach 
nes and Electromechanics, n.º 3, 1978. 

7-—N. HANCOCK — «Matrix Analysis of Electric Machi- 
nery»; Pergamon Press. 

8 — T. A. LIPO and P. €C. KRAUSE — «Stability Analysjs of 
a Rectifier Inverter Induction Motor Drive»; 1. E. E. E. 
Transactions on Power Apparatus and Systems, vol. 
PAS-88, n.º 1, January 1969. 

9-—M. S. GARRIDO and E. GUDEFIN — «ÉEquations des 
circujts à commutation non linéiques»; Revue E, vol VII, 
nº 3, 1972. 


ACKNOWLEDGEMENT 


The authors are grateful to Professor M. S. Garrido 
for the profitable discussion on the subject of this 
paper. 


101 


STRIP GRATING LEAKY-WAVE ANTENNA 


A.M. BARBUSA e J. FIGANIER 


Departamento de Engenharia Eelectrotécnica e de Computadores, Secção de Propagação 
e Radiação e Centro de Análise e Processamento de Sinais, IST 


ABSTRACT 


The use as a leaky-wave antenna of a periodie metallic 
strip grating over a dielectric slab backed by a perfectly 
conducting plane is discussed. Typical seis of parameters 
are given, 


RESUMO 


Discute-se a aplicação como antena “leaky-wave” de 
uma estrutura constituída por uma placa dieléctrica sobre 
um plano metálico condutor e coberta por uma grelha de 
tiras metálicas. Apresentam-se valores típicos dos parâme- 
tros desta estrutura, 


INTRODUCTION 


A grounded dielectric slab covered by a strip grating has 
been studied earlier by Sigelmann and Ishimaru!, Sigel- 
mann?, Jacobsen? and Balling?. In their analysis they have 
used approximations which are valid for narrow slots!-? or 
narrow strips**. Sigelmann discussed only solutions of the 
modal equation of the surface wave type. Jacobsen discus- 
sed also the leaky wave solutions and gave numerical re- 
sults for structures where the dielectric slab has a high die- 
lectric constant. Hence the fundamental spatial harmonic is 
a slow wave and radiation is expected to be due only to 
higher order harmonies. Balling dealt with a structure of 
the same type and examined the range of parameters for a 
viable antenna with narrow slots radiating by an higher 
order harmonic, 

More recently several leaky-wave antennas have been 
proposed based on different types of dielectric waveguides 
perturbed by strip gratings*8. In these structures the domi- 
nant mode is a slow wave and radiation is due to higher 
order harmonics, namely n = | 1. When the strips are 
narrow and the perturbation due to them is small the gra- 
ting effect may be analysed in terms of a linear array. An 
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analysis based on a spectral domain approach has been 
presented for a leaky wave antenna composed of an array 
of thin metallic strips on a rectangular dielectric rod”. 

Elsewhere!O the authors have studied the wave guiding 
characteristics of a metallc strip grating running parallel to 
an impedance plane using a different method which is ap- 
plicable to a structure with an arbitrary value of the slot 
width-period ratio and have shown that this structure sup- 
ports leaky waves. However the case considered refers to a 
structure with homogeneous dielectric which is not ade- 
quate for a practical antenna. 

In this note a structure where the metallic strips are pla- 
ced periodically along a grounded dielectric slab is analy- 
sed. As was to be expected this structure also supports 
leaky waves and if the parameters are suitably chosen the 
design of a practical leaky-wave antenna may be based on 
It. 


ANALYSIS AND MODAL EQUATION 


A time dependence exp (ju t) is assumed and field pro- 
pagation along the z axis (Fig. 1) is represented by the 
tactor exp (-J À 2). Only TM propagation with respect to 
the z axis is studied, 

In mathematical terms we are dealing with a boundary 
value problem with mixed type boundary conditions which 
is formulated as a modified Riemann-Hilbert problem!º. 
The details of the calculations are omitted since we follow 
closely Ref. 10. Due to the presence of two dielectrics we 
have Bo two sets of transverse wavenumbers ha, sa 
= (Ma - Ma e, ko = 0 (615 no)? and Aq = 
=A-n7f, 2 being the period of the structure. To 
simplify the notation we use henceforth h, and h respectively 
for ho and ho,. Itis assumed that Im(hn,) < o,forn > 1. 
This condition is sufficient to ensure the convergence of the 
series for the potential but allows the existence of leaky 


wave solutions!º. 


Original recebido para publicação em 30/11/84. 


The calculations lead to the following expression for the 
modal equation 


f(h) = 


=fe[1 + exp(-2jhy b)J + Pl — exp( 2jh b)]) 601 — 1) 


o Am (h) 
E A sin(m 81) 
m = 
, + sm(môB,) 
e a ne (1) 
n m m 


where the prime stands for omission of the terms n = o and 
m=0;0,=md/L,e=€e/€ and 


Gm =) Ed (1 +ed)hy [exp(-2jh;b) = 1] 
[P, (cos 81) - P (cos8,)] 


In the last expression, P | is the Legendre polynomial of 
order m. The functions u, are the solutions of an infinite 
system of algebraic equations whose matrixis A = [A |]. 
This system is the same as for the case of a single dielectric 
structure!O except for the expression of the perturbation term 
which now depends on e,. 

Solutions of equation (1) in the upper half of the h — plane 
correspond to improper modes (leaky wave modes) and 
solutions in the lower half of the h — plane correspond to 
proper modes; those on the negative imaginary axis are 
surface wave modes. 


Fig. | 


Geometry of the structure 


NUMERICAL RESULTS AND CONCLUSIONS 


The modal equation was solved numerically using the 
method mentioned in Ref. 10. Numerical results have been 
obtained for structures with periods 24 < AG /2 Ag 
being the free space wavelength. Iíthe propagation constant 
of the fundamental harmonic is near to ks the condition 
24 < Ag) 2ensures that the propagation constants of the 
higher order harmonics lie outside the visible region. Thus 
a single beam radiating almost endfire may be obtained. 

The numerical results show that only structures with small 
slot widths, in the order of 0.054, support a dominant mode 
with leakage constant not too high. This is an expected 
result taking into account that the strip currents are per- 
pendicular to the edges. Larger slot widths give rise to atte- 
nuation constants too large to be appropriate for an an- 
tenna with a narrow beam and a good sidelobe level. 

In Figs. 2 and 3 the longitudinal phase velocity (v) and 
the longitudinal attenuation constant (a) of the dominant 
leaky-wave mode are given for some sets of parameters. 
Curve | applies to a single dielectric (air) structure and is 
given for comparison. Fig. 2a shows that a dominant mode 
with v/c = lis possible only with a small e, (curves 2 and 


3). Increasing the relative dielectric constant (curve 4) gives 
rise to slow modes and radiation is no longer possible due 
to the fundamental harmonic. Moreover the increase in e, 
reduces the useful bandwidth of the structure. In Figures 2b 
and 3b, curves 2 and 3 differ only by the slot width, (curve | 
is the same as curve 2 in Figures 2a and 3a and is repeated 
for comparison), showing that the leakage may be control- 
led almost without affecting the phase velocity. 


Fig. 2a Phase velocity of the first leaky mode 


£—d = 0,05 
b=4 
curve | 2 3 4 
é; 1.0 1.05 1.1 1.3 


Fig. 2b Phase velocity of the first leaky mode 


6. = 1.05 
curve | 2 3 
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Fig. 3a Attenuation constant of the first leaky mode. 
Same parameters as in Fig. 2a. 
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Fig. 3b Attenuation constant of the first leaky mode. 
Same parameters as in Fig. 2b. 
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Fig. 4 Phase velocity and attenuation constant of the 
second leaky mode £ —- d = 00254, b = 2.5€, 
e =3 


The thickness b of the dielectric also controlls the radia- 
tion properties of the structure; for a small b the first leaky 
mode is dominant; for a greater value of b the leakage 
constant of the second leaky mode decreases (Fig. 3b) and 
may become comparable to the first leaky mode constant. 
Thus single mode operation may become impossible. For 
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structures with small slot widths we may estimate the al- 
lowed values of b for single mode operation by taking into 
account that the second leaky mode tends to the TM, mode 
of the parallel plate waveguide with the same dielectric and 
thickness. 

The operation of the strip grating in the second leaky mode 
may also be of some interest. For a structure with a higher 
dielectric constant, say e, = 3, keeping the same frequency 
and geometrical parameters, the first mode is no longer a 
leaky mode but becomes an improper surface mode * (with 
root locus in the positive imaginary axis of the h-plane) whilst 
the second leaky mode becomes a fast wave with phase 
velocity (Fig. 4) suitable for frequency scanning operation 
although only for a smaller frequency range. 

The radiation of a structure of this class in the presence 
of a suitable primary source has been studied and the re- 
sults, namely the radiation patterns, are reported elsew- 
here!2. 
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This mode does not contribute to the radiation field in the 
steepest descent representation 


TOPOGRAF IA 


DOMINGOS MOURA 


TÉCNICAS 
DE ; 
ALTA TENSÃO 


8" .. ] Ra Ea o CN 
= 


TABELAS 
TÉCNICAS 


É 


TO ABECASIS MANZANARES 


=. 
.... 


AIJ) 


LST DAS 
e am 


NE 
o. a 


- ESCOAMENTOS LÍQUIDOS 


a... 
o... 


a... .. 
ad erre 


LE: 
rea 


TÉCNICA *- astst » TREÉ 


Pós ddddo 
SICOROCANAANANNTAE: 
NOEDANANaAaaa 


oa... 
NESONERR. 

ALTAS A 

inosesases 
an. é 
a... 
a... 

ei... ia 

a EEE 


REVISTA DE ENGENHARIA DA ASSOCIAÇÃO DOS ESTUDANTES 
DO INSTITUTO SUPERIOR TECNICO - LISBOA 


Acerca do artigo , 
SOBRE A REPRESENTAÇÃO ANALÍTICA 
DE RESULTADOS EXPERIMENTAIS 


“TÉCNICA” 42:23-26, n.º 460, Out. 1980 


P. DE VARENNES E MENDONÇA 


O problema da regressão linear clássica põe-se no espaço 
cartesiano n-dimensional e visa determinar a . quação do 


hiperplano 
Y=ajtayX)+a;X3;+...+a, Xn 
a partir da amostra (y;; Xo. X3jo co. Xni). 


comi=1,2.3..m > n. 


Usa-se de ordinário o método dos mínimos quadrados, 
procedimento justificado por os seus estimadores dos pará- 
metros a, serem BLUE, satisfeitas que se encontrem as con- 
dições do teorema de Gauss-Markov [v. g.:3, p. 861. 

Para valores altos de n o cálculo numérico é assaz tra- 
balhoso, facto que leva na habitual prática hodierna a utili- 
zar computadores. 

Não assim para n = 4. Se muitas calculadoras electró- 
nicas de bolso resolvem com simplicidade o problema para 
n = 2, também algumas há com as quais as constantes a, 
a, c ay da equação do plano de regressão (n = 3) 


1) Y=aj+a X,+ay X; 
e até mesmo asa|,,a,,a,e as de 
2) Y=aj+aX,+ayX, + ay Xs 
encontram fácil determinação. 

Quanto a 2) pode citar-se a HP - 41C equipada com 
a ROM dita módulo STAT I[2,p. 38]: (quanto a 1). as 
TI-58, TI-58€ e TI-59 com o módulo n. 2. 


No caso do exemplo apresentado no artigo em foco pre- 
tende-se achar a equação da parábola, 


é 
3) Y=a+a X,+ as X.. 
projecção da intersecção do plano |) com a superficie cilin- 


É E E: 
drica X, = X5. 


Havendo razões para seleccionar as variáveis explicativas, 
com consequente redução do valor de n, adoptar-se-á um 
método racional, como seja o algoritmo STRAP [1]. A 
hipótese (cuja admissão enferma sempre de subjectividade) 


Original recebido para publicação em 29/4/81. 
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exigida pelo processo proposto no artigo de existir um ponto 
cujas coordenadas são conhecidas exactamente e pelo qual 
su à passar 3) é todavia em geral inconveniente, por- 
quanto diminui arbitrária e simultaneamente o número de 
graus de liberdade e o coeficiente de determinação. 

Com uma das calculadoras TI citadas, o processo cor- 
recto é o seguinte [4,p. 3-8e5-14]: PGMOS E' x, A x) By, 
Cx Axo, By; Cx Axo Bym C PGM 
IS A-aB- a C- as. 

Para os dados do exemplo obtém-se a, = 0.998 1964286, 
a, = — 0.0019023214 e a; = 0.0000272321. 


Querendo comparar os valores estimados pela parábola 
de regressão com os dados, continua-se assim: (FIX 4) 
OA x? B'— 0.99822 A' x2 B' — 0.99454 A' x2 B' — 
—- 0.9910 6 A' x2 B'— 0.9878 8 A' x? B' — 0,9847 
IO A' x? B' — 0.9819. O resultado denota um ajustamento 
perfeito e torna dispensável qualquer outro método de ava- 
liação da sua qualidade. Não há de resto dificuldade al- 
guma em calcular o coeficiente de determinação, já que 
basta premir a tecla D para se obter R? = 1.0000, 

É essa perfeição que mascara a inconveniência da intro- 
dução da mencionada hipótese, inconveniência que qualquer 
exemplo onde Rº se afaste um tanto da unidade logo evi- 
dencia. 
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ERRATA 
pag. coluna linha onde se lê deve ler-se 
55 2 9 ««. 2: à contribuições do ... «.. duas contribuições: a do ... 
56 1 52 sos 108 CF sxs sas 08º 
RMR - 10 
59 1 3 o. By ES... cad 4 
29 l 19-20 Falta uma linha -.-.« das descontinuidades estão re 
lacionadas não sô com ... 

58 Quadro C - 2a coluna > 56 MPa > 56 000 MPa 

Quadro D - 2a coluna 250 MPa > 250 MPa 

Quadro D - 2a coluna 65º >: 65” 

Quadro E - 6a coluna < 5 mm >” 5 mm 
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